diff --git a/DESCRIPTION b/DESCRIPTION index d5332d5..aceed93 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ResistanceGA Title: Optimize resistance surfaces using genetic algorithms -Version: 3.1-3 +Version: 4.0 Author: Bill Peterman Maintainer: Bill Peterman Description: This package contains functions to optimize @@ -14,6 +14,7 @@ Depends: raster Imports: ggplot2, + ggExtra, akima, plyr, dplyr, @@ -23,7 +24,7 @@ Imports: lme4 (>= 1.0), Matrix, MuMIn, - smoothie + spatstat Suggests: knitr, rmarkdown diff --git a/NAMESPACE b/NAMESPACE index c9f7e6f..21924b9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -20,9 +20,12 @@ export(Run_CS) export(Run_gdistance) export(SS_optim) export(SS_optim.scale) +export(To.From.ID) +export(all_comb) export(gdist.prep) export(k.smooth) export(lower) +export(mlpe_rga) import(GA) import(gdistance) import(ggplot2) @@ -37,6 +40,8 @@ importFrom(dplyr,left_join) importFrom(dplyr,mutate) importFrom(dplyr,summarise) importFrom(dplyr,tally) +importFrom(ggExtra,ggMarginal) +importFrom(ggExtra,removeGrid) importFrom(grDevices,dev.off) importFrom(grDevices,tiff) importFrom(grDevices,topo.colors) @@ -49,8 +54,11 @@ importFrom(plyr,create_progress_bar) importFrom(plyr,ldply) importFrom(plyr,progress_text) importFrom(plyr,rbind.fill) -importFrom(smoothie,kernel2dsmooth) +importFrom(spatstat,as.im) +importFrom(spatstat,as.matrix.im) +importFrom(spatstat,blur) importFrom(stats,AIC) +importFrom(stats,as.formula) importFrom(stats,lm) importFrom(stats,logLik) importFrom(stats,qqline) @@ -58,7 +66,10 @@ importFrom(stats,qqnorm) importFrom(stats,resid) importFrom(stats,residuals) importFrom(stats,runif) +importFrom(stats,sigma) +importFrom(utils,combn) importFrom(utils,file_test) +importFrom(utils,menu) importFrom(utils,read.csv) importFrom(utils,read.delim) importFrom(utils,read.table) diff --git a/NEWS b/NEWS index ce93652..e0097c6 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,13 @@ +ResistanceGA 4.0-0 +-------- +* This is the official release of the development version that includes the `all_comb` function, as well as a generic `mlpe_rga` function for flexibly fitting MLPE mixed effects models. + +ResistanceGA 3.1-4 +-------- +* Updated Plot.trans function to generate marginal plots for response curve transformations. This uses the ggExtra package +* Have implemented an `all_combs` function to facilitate the running of a full analysis from a single function: replicate runs of single and multisurface combinations, bootstrap resampling of results. +* Gaussian smoothing now done with spatstat rather than smoothie package + ResistanceGA 3.1-3 -------- * Default setting for k.value in GA.prep is now 2, meaning that AICc is calculated based on the number of parameters optimized, plus 1 for the mixed effects model intercept. diff --git a/R/Bootstrap_function.R b/R/Bootstrap_function.R index fa68a9b..2799280 100644 --- a/R/Bootstrap_function.R +++ b/R/Bootstrap_function.R @@ -8,14 +8,21 @@ #' @param sample.prop Proportion of observations to be sampled each iteration (Default = 0.75) #' @param iters Number of bootstrap iterations to be conducted #' @param obs Total number of observations (populations or individuals) in your original analysis +#' @param rank.method What metric should be used to rank models during bootstrap analysis? c('AIC', 'AICc', 'R2', 'LL'). Default = 'AICc' #' @param genetic.mat Genetic distance matrix without row or column names. -#' @return A data frame reporting the average model weight, average rank, number of time a model was the top model in teh set, and the frequency a model was best. +#' @return A data frame reporting the average model weight, average rank, number of time a model was the top model in the set, and the frequency a model was best. +#' +#' @details This is a 'pseudo'bootstrap procedure that subsamples distance and genetic matrices, refits the MLPE model for each surface. AICc is calculated based on the number of parameters specified. Ranking of models during the bootstrap analysis is based on the specified \code{rank.method}, which defaults to 'AICc'. The objective of this procedure is to identify the surfaces that is top ranked across all bootstrap iterations. #' -#' @details This is a 'pseudo'bootstrap procedure that subsamples distance and genetic matrices, refits the MLPE model for each surface, and assesses the model AICc. AICc is calculated based on the number of parameters specified. - #' @export #' @author Bill Peterman -#' @usage Resist.boot(mod.names, dist.mat, n.parameters, sample.prop, iters, obs, genetic.mat) +#' @usage Resist.boot(mod.names, +#' dist.mat, +#' n.parameters, +#' sample.prop, +#' iters, obs, +#' rank.method, +#' genetic.mat) #' Resist.boot <- @@ -25,6 +32,7 @@ Resist.boot <- sample.prop = 0.75, iters, obs, + rank.method = 'AICc', genetic.mat) { options(warn = -1) @@ -62,7 +70,7 @@ Resist.boot <- ) mod.aic <- data.frame(mod.names[j], n.parameters[j], AICc) - names(mod.aic) <- c("surface", "k", "AICc", "R2m") + names(mod.aic) <- c("surface", "k", "AIC","AICc", "R2m", "LL") AICc.tab[[j]] <- mod.aic progress_bar$step() @@ -71,12 +79,24 @@ Resist.boot <- # Calculate Delta AICc, weight, and rank AICc.tab <- plyr::ldply(AICc.tab, "identity") - AICc.tab <- AICc.tab %>% plyr::mutate(., delta = AICc - min(AICc)) %>% + AICc.tab <- AICc.tab %>% dplyr::mutate(., AIC = AIC) %>% + dplyr::mutate(., AICc = AICc) %>% + plyr::mutate(., delta = AICc - min(AICc)) %>% dplyr::mutate(., weight = (exp(-0.5 * delta)) / sum(exp(-0.5 * delta))) %>% - dplyr::mutate(., rank = rank(AICc)) %>% dplyr::mutate(., iteration = i) %>% + dplyr::mutate(., LL = LL) %>% dplyr::mutate(., R2m = R2m) + if(rank.method == 'LL') { + AICc.tab <- dplyr::mutate(AICc.tab, rank = dense_rank(desc(LL))) + } else if (rank.method == 'AIC') { + AICc.tab <- dplyr::mutate(AICc.tab, rank = dense_rank(AIC)) + } else if(rank.method == 'R2'){ + AICc.tab <- dplyr::mutate(AICc.tab, rank = dense_rank(desc(R2m))) + } else { + AICc.tab <- dplyr::mutate(AICc.tab, rank = dense_rank(AICc)) + } + AIC.tab.list[[i]] <- AICc.tab } # Close iteration loop @@ -85,9 +105,12 @@ Resist.boot <- boot.avg <- group.list %>% dplyr::summarise(., - avg.R2m = mean(R2m), + avg.AIC = mean(AIC), + avg.AICc = mean(AICc), avg.weight = mean(weight), - avg.rank = mean(rank)) %>% + avg.rank = mean(rank), + avg.R2m = mean(R2m), + avg.LL = mean(LL)) %>% arrange(., avg.rank) Freq_Percent <- @@ -121,7 +144,7 @@ To.From.ID <- function(POPS) { colnames(tmp2) <- c("pop1", "pop2") as.numeric(tmp2$pop1) as.numeric(tmp2$pop2) - ID <- arrange(tmp2, as.numeric(pop1), as.numeric(pop2)) + ID <- plyr::arrange(tmp2, as.numeric(pop1), as.numeric(pop2)) # ID<-tmp2[with(tmp2, order(pop1, pop2)), ] p1 <- ID[POPS - 1, 1] p2 <- ID[POPS - 1, 2] @@ -156,7 +179,11 @@ boot.AICc <- function(response, resistance, ID, ZZ, k) { fit.mod <- mkMerMod(environment(dfun), opt, mod$reTrms, fr = mod$fr) R.sq <- MuMIn::r.squaredGLMM(fit.mod)[[1]] mod.AIC <- AIC(fit.mod) + LL <- logLik(fit.mod) AICc <- mod.AIC + ((2 * k * (k + 1)) / (length(response) - k - 1)) - out <- data.frame(AICc = AICc, R2m = R.sq) + out <- data.frame(AIC = mod.AIC, + AICc = AICc, + R2m = R.sq, + LL = LL) return(out) } \ No newline at end of file diff --git a/R/Combine_Surfaces.R b/R/Combine_Surfaces.R index 487414d..6665875 100644 --- a/R/Combine_Surfaces.R +++ b/R/Combine_Surfaces.R @@ -25,7 +25,14 @@ #' #' The Distance transformation sets all values equal to one. Because of the flexibility of the Ricker function to take a monomolecular shape (try \code{Plot.trans(PARM=c(10,100), Resistance=c(1,10), transformation="Ricker")} to see this), whenever a shape parameter >6 is selected in combination with a Ricker family transformation, the transformation reverts to a Distance transformation. In general, it seems that using a combination of intermediate Ricker and Monomolecular transformations provides the best, most flexible coverage of parameter space. #' @return R raster object that is the sum all transformed and/or reclassified resistance surfaces provided -#' @usage Combine_Surfaces(PARM, CS.inputs, gdist.inputs, GA.inputs, out, File.name, rescale, p.contribution) +#' @usage Combine_Surfaces(PARM, +#' CS.inputs, +#' gdist.inputs, +#' GA.inputs, +#' out, +#' File.name, +#' rescale, +#' p.contribution) #' @export #' @author Bill Peterman Combine_Surfaces <- diff --git a/R/Combine_Surfaces_scale.R b/R/Combine_Surfaces_scale.R index a7e5303..e21eaf5 100644 --- a/R/Combine_Surfaces_scale.R +++ b/R/Combine_Surfaces_scale.R @@ -25,7 +25,14 @@ #' #' The Distance transformation sets all values equal to one. Because of the flexibility of the Ricker function to take a monomolecular shape (try \code{Plot.trans(PARM=c(10,100), Resistance=c(1,10), transformation="Ricker")} to see this), whenever a shape parameter >6 is selected in combination with a Ricker family transformation, the transformation reverts to a Distance transformation. In general, it seems that using a combination of intermediate Ricker and Monomolecular transformations provides the best, most flexible coverage of parameter space. #' @return R raster object that is the sum all transformed and/or reclassified resistance surfaces provided -#' @usage Combine_Surfaces.scale(PARM, CS.inputs, gdist.inputs, GA.inputs, out, File.name, rescale, p.contribution) +#' @usage Combine_Surfaces.scale(PARM, +#' CS.inputs, +#' gdist.inputs, +#' GA.inputs, +#' out, +#' File.name, +#' rescale, +#' p.contribution) #' @export #' @author Bill Peterman Combine_Surfaces.scale <- @@ -150,6 +157,9 @@ Combine_Surfaces.scale <- parm <- PARM[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + 1])] + if(parm[4] < 0.5) { + parm[4] <- 0.000123456543210 + } rast <- k.smooth(raster = r[[i]], sigma = parm[4], diff --git a/R/GA_prep.R b/R/GA_prep.R index 76ae7eb..393b623 100644 --- a/R/GA_prep.R +++ b/R/GA_prep.R @@ -3,13 +3,13 @@ #' This function prepares and compiles objects and commands for optimization with the GA package #' #' @param ASCII.dir Directory containing all raster objects to optimized. If optimizing using least cost paths, a RasterStack or RasterLayer object can be specified. -#' @param Results.dir If a RasterStack is provided in place of a directory containing .asc files for ASCII.dir, then a directory to export optimization results must be specified. It is critical that there are NO SPACES in the directory, as this will cause the function to fail. +#' @param Results.dir If a RasterStack is provided in place of a directory containing .asc files for ASCII.dir, then a directory to export optimization results must be specified. It is critical that there are NO SPACES in the directory, as this will cause the function to fail. If using the \code{\link[ResistanceGA]{all_comb}} function, specify \code{Results.dir} as "all.comb". #' @param min.cat The minimum value to be assessed during optimization of of categorical resistance surfaces (Default = 1e-04) #' @param max.cat The maximum value to be assessed during optimization of of categorical resistance surfaces (Default = 2500) #' @param max.cont The maximum value to be assessed during optimization of of continuous resistance surfaces (Default = 2500) -#' @param min.scale The minimum scaling parameter value to be assessed during optimization of resistance surfaces with kernel smoothing (Default = 1) -#' @param max.scale The maximum scaling parameter value to be assessed during optimization of resistance surfaces with kernel smoothing (Default = 0.25 * nrows in the raster surface) -#' @param cont.shape A vector of hypothesized relationships that each continuous resistance surface will have in relation to the genetic distance reposnse (Default = NULL; see details) +#' @param min.scale The minimum scaling parameter value to be assessed during optimization of resistance surfaces with kernel smoothing (Default = 0.01). See details +#' @param max.scale The maximum scaling parameter value to be assessed during optimization of resistance surfaces with kernel smoothing (Default = 0.1 * maximum dimension of the raster surface) +#' @param cont.shape A vector of hypothesized relationships that each continuous resistance surface will have in relation to the genetic distance response (Default = NULL; see details) #' @param select.trans Option to specify which transformations are applied to continuous surfaces. Must be provided as a list. "A" = All, "M" = Monomolecular only, "R" = Ricker only. See Details. #' @param method Objective function to be optimized. Select "AIC", "R2", or "LL" to optimize resistance surfaces based on AIC, variance explained (R2), or log-likelihood. (Default = "LL") #' @param scale Logical. To optimize a kernel smoothing scaling parameter during optimization, set to TRUE (Default = FALSE). See Details below. @@ -38,7 +38,7 @@ #' @param run Number of consecutive generations without any improvement in AICc before the GA is stopped (Default = 25) #' @param keepBest A logical argument specifying if best solutions at each iteration should be saved (Default = TRUE) #' @param seed Integer random number seed to replicate \code{ga} optimization -#' @param quiet Logical. If TRUE, AICc and step run time will not be printed to the screen after each step. Only \code{ga} summary information will be printed following each iteration. Default = FALSE +#' @param quiet Logical. If TRUE, the objective function and step run time will not be printed to the screen after each step. Only \code{ga} summary information will be printed following each iteration. (Default = FALSE) #' @return An R object that is a required input into optimization functions #' #' @details Only files that you wish to optimize, either in isolation or simultaneously, should be included in the specified \code{ASCII.dir}. If you wish to optimize different combinations of surfaces, different directories contaiing these surfaces must be created. @@ -47,6 +47,8 @@ #' #' \code{scale.surfaces} can be used to specify which surfaces to apply kernel smoothing to during multisurface optimization. For example, \code{scale.surfaces = c(1, 0, 1)} will result in the first and third surfaces being optimized with a kernel smoothing function, while the second surface will not be scaled. The order of surfaces will match either the order of the raster stack, or alphabetical order when reading in from a directory. #' +#' \code{min.scale} defaults to a minimum of 0.01. During optimization, whenever the scaling factor (sigma) is less than 0.5, ResistanceGA will not apply scaling. In this way, it is possible for a surface to not be scaled. +#' #' The Default for \code{k.value} is 2, which sets k equal to the number of parameters optimized, plus 1 for the intercept term. Prior to version 3.0-0, \code{k.value} could not be specified by the user and followed setting 2, such that k was equal to the number of parameters optimized plus the intercept term. #' #' \code{cont.shape} can take values of "Increase", "Decrease", or "Peaked". If you believe a resistance surface is related to your response in a particular way, specifying this here may decrease the time to optimization. \code{cont.shape} is used to generate an initial set of parameter values to test during optimization. If specified, a greater proportion of the starting values will include your believed relatiosnship. If unspecified (the Default), a completely random set of starting values will be generated. @@ -122,7 +124,7 @@ GA.prep <- function(ASCII.dir, scale <- NULL } - if (!is.null(Results.dir)) { + if (!is.null(Results.dir) & (Results.dir != 'all.comb')) { TEST.dir <- !file_test("-d", Results.dir) if (TEST.dir == TRUE) { stop("The specified 'Results.dir' does not exist") @@ -171,18 +173,26 @@ GA.prep <- function(ASCII.dir, stop("The 'scale.surfaces' vector is not the same length as the number of layers") } - if ("Results" %in% dir(Results.dir) == FALSE) - dir.create(file.path(Results.dir, "Results")) - # dir.create(file.path(ASCII.dir, "Results"),showWarnings = FALSE) - Results.DIR <- paste0(Results.dir, "Results/") - if ("tmp" %in% dir(Results.dir) == FALSE) - dir.create(file.path(Results.dir, "tmp")) - # dir.create(file.path(Results.dir, "tmp"),showWarnings = FALSE) - Write.dir <- paste0(Results.dir, "tmp/") - if ("Plots" %in% dir(Results.DIR) == FALSE) - dir.create(file.path(Results.DIR, "Plots")) - # dir.create(file.path(Results.dir, "tmp"),showWarnings = FALSE) - Plots.dir <- paste0(Results.DIR, "Plots/") + if(Results.dir != 'all.comb') { + if ("Results" %in% dir(Results.dir) == FALSE) + dir.create(file.path(Results.dir, "Results")) + Results.DIR <- paste0(Results.dir, "Results/") + + if ("tmp" %in% dir(Results.dir) == FALSE) + dir.create(file.path(Results.dir, "tmp")) + Write.dir <- paste0(Results.dir, "tmp/") + + if ("Plots" %in% dir(Results.DIR) == FALSE) + dir.create(file.path(Results.DIR, "Plots")) + Plots.dir <- paste0(Results.DIR, "Plots/") + } + + if(Results.dir == 'all.comb') { + Results.DIR <- NULL + Write.dir <- NULL + Plots.dir <- NULL + } + # Determine total number of parameters and types of surfaces included parm.type <- data.frame() @@ -202,17 +212,19 @@ GA.prep <- function(ASCII.dir, ) } - if (n.levels == 2 & !is.null(scale) & scale.surfaces[i]==1) { + if (n.levels == 2 & + !is.null(scale) & + scale.surfaces[i]==1) { parm.type[i, 1] <- "cont" parm.type[i, 2] <- 4 parm.type[i, 3] <- names[i] if (is.null(min.scale)) { - min.scale <- 0.3 - + min.scale <- 0.01 } + if (is.null(max.scale)) { - max.scale <- max(dim(r[[i]])) / 4 + max.scale <- floor(max(dim(r[[i]])) * 0.1) } min.list[[i]] <- @@ -246,23 +258,23 @@ GA.prep <- function(ASCII.dir, parm.type[i, 2] <- 4 if (is.null(min.scale)) { - min.scale <- 0.3 - + min.scale <- 0.01 + } if (is.null(max.scale)) { - max.scale <- max(dim(r[[i]])) + max.scale <- floor(max(dim(r[[i]])) * 0.1) } min.list[[i]] <- - c(1, 0.001, 0.001, min.scale) # eq, shape/gaus.opt, max, gaus.sd + c(1, 0.5, 0.001, min.scale) # eq, shape/gaus.opt, max, gaus.sd max.list[[i]] <- c(9.99, 15, max.cont, max.scale) } else { parm.type[i, 2] <- 3 min.list[[i]] <- - c(1, 0.001, 0.001) # eq, shape/gaus.opt, max - - max.list[[i]] <- c(9.99, 15, max.cont) + c(1, 0.5, 0.001) # eq, shape/gaus.opt, max + + max.list[[i]] <- c(9.99, 14.5, max.cont) } parm.type[i, 3] <- names[i] @@ -298,6 +310,7 @@ GA.prep <- function(ASCII.dir, pop.size = pop.size, max = max.cont, scale = scale, + min.scale = min.scale, max.scale = max.scale ) @@ -316,6 +329,7 @@ GA.prep <- function(ASCII.dir, pop.size = pop.size, max = max.cont, scale = scale, + min.scale = min.scale, max.scale = max.scale ) cont.shape <- cont.shape#[-1] @@ -332,6 +346,7 @@ GA.prep <- function(ASCII.dir, pop.size = pop.size, max = max.cont, scale = scale, + min.scale = min.scale, max.scale = max.scale ) } else { @@ -349,43 +364,117 @@ GA.prep <- function(ASCII.dir, Min.Max <- 'max' } - list( - parm.index = parm.index, - ga.min = ga.min, - ga.max = ga.max, - select.trans = eqs, - scale = scale, - scale.surfaces = scale.surfaces, - surface.type = surface.type, - parm.type = parm.type, - Resistance.stack = r, - n.layers = n.layers, - layer.names = names, - pop.size = pop.size, - min.list = min.list, - max.list = max.list, - SUGGESTS = SUGGESTS, - ASCII.dir = ASCII.dir, - Results.dir = Results.DIR, - Write.dir = Write.dir, - Plots.dir = Plots.dir, - type = type, - pcrossover = pcrossover, - pmutation = pmutation, - crossover = crossover, - maxiter = maxiter, - run = run, - keepBest = keepBest, - population = population, - selection = selection, - mutation = mutation, - parallel = parallel, - pop.mult = pop.mult, - percent.elite = percent.elite, - Min.Max = Min.Max, - method = method, - k.value = k.value, - seed = seed, - quiet = quiet - ) + if(Results.dir != "all.comb") { + list( + parm.index = parm.index, + ga.min = ga.min, + ga.max = ga.max, + select.trans = eqs, + scale = scale, + scale.surfaces = scale.surfaces, + surface.type = surface.type, + parm.type = parm.type, + Resistance.stack = r, + n.layers = n.layers, + layer.names = names, + pop.size = pop.size, + min.list = min.list, + max.list = max.list, + SUGGESTS = SUGGESTS, + ASCII.dir = ASCII.dir, + Results.dir = Results.DIR, + Write.dir = Write.dir, + Plots.dir = Plots.dir, + type = type, + pcrossover = pcrossover, + pmutation = pmutation, + crossover = crossover, + maxiter = maxiter, + run = run, + keepBest = keepBest, + population = population, + selection = selection, + mutation = mutation, + parallel = parallel, + pop.mult = pop.mult, + percent.elite = percent.elite, + Min.Max = Min.Max, + method = method, + k.value = k.value, + seed = seed, + quiet = quiet + ) + } else { + list( + parm.index = parm.index, + ga.min = ga.min, + ga.max = ga.max, + select.trans = eqs, + scale = scale, + scale.surfaces = scale.surfaces, + surface.type = surface.type, + parm.type = parm.type, + Resistance.stack = r, + n.layers = n.layers, + layer.names = names, + pop.size = pop.size, + min.list = min.list, + max.list = max.list, + SUGGESTS = SUGGESTS, + ASCII.dir = ASCII.dir, + Results.dir = Results.DIR, + Write.dir = Write.dir, + Plots.dir = Plots.dir, + type = type, + pcrossover = pcrossover, + pmutation = pmutation, + crossover = crossover, + maxiter = maxiter, + run = run, + keepBest = keepBest, + population = population, + selection = selection, + mutation = mutation, + parallel = parallel, + pop.mult = pop.mult, + percent.elite = percent.elite, + Min.Max = Min.Max, + method = method, + k.value = k.value, + seed = seed, + quiet = quiet, + inputs = list( + ASCII.dir = ASCII.dir, + Results.dir = Results.dir, + min.cat = min.cat, + max.cat = max.cat, + max.cont = max.cont, + min.scale = min.scale, + max.scale = max.scale, + cont.shape = cont.shape, + select.trans = select.trans, + method = method, + scale = scale, + scale.surfaces = scale.surfaces, + k.value = k.value, + pop.mult = pop.mult, + percent.elite = percent.elite, + type = type, + pcrossover = pcrossover, + pmutation = pmutation, + maxiter = maxiter, + run = run, + keepBest = keepBest, + population = population, + selection = selection, + crossover = crossover, + mutation = mutation, + pop.size = pop.size, + parallel = parallel, + seed = seed, + quiet = quiet + ) + ) + } + } \ No newline at end of file diff --git a/R/MLPE.R b/R/MLPE.R index 5c5f850..d90c2b2 100644 --- a/R/MLPE.R +++ b/R/MLPE.R @@ -1,7 +1,7 @@ -# Run Mixed effects models, recovery parameter estimates +# Run Mixed effects models, recover parameter estimates #' Run maximum likelihood population effects mixed effects model (MLPE) #' -#' Runs MLPE as detailed by Clarke et al. (2002). This function will run the model and return glmer object +#' Runs MLPE as detailed by Clarke et al. (2002). This function will run the model and return lmer object #' #' @param resistance Path to pairwise resistance distance matrix (resistances.out) from CS results. Alternatively, provide the pairwise resistances created from optimizing with `gdistance` (result of Run_gdistance). #' @param pairwise.genetic Lower half of pairwise genetic distance matrix @@ -14,7 +14,12 @@ #' @export #' @author Bill Peterman -#' @usage MLPE.lmm(resistance, pairwise.genetic, REML, ID, ZZ) +#' @usage MLPE.lmm(resistance, +#' pairwise.genetic, +#' REML = FALSE, +#' ID = NULL, +#' ZZ = NULL, +#' scale = TRUE) #' @references Clarke, R. T., P. Rothery, and A. F. Raybould. 2002. Confidence limits for regression relationships between distance matrices: Estimating gene flow with distance. Journal of Agricultural, Biological, and Environmental Statistics 7:361-372. MLPE.lmm <- diff --git a/R/MS_optim.R b/R/MS_optim.R index 56743e0..6a666a7 100644 --- a/R/MS_optim.R +++ b/R/MS_optim.R @@ -205,11 +205,34 @@ MS_optim <- function(CS.inputs = NULL, cd.list <- list(as.matrix(cd)) names(cd.list) <- NAME + AICc.tab <- data.frame(surface = NAME, + obj = multi.GA_nG@fitnessValue, + k = k, + AIC = aic, + AICc = AICc, + R2m = fit.stats[[1]], + R2c = fit.stats[[2]], + LL = LL) + + colnames(AICc.tab) <- + c( + "Surface", + paste0("obj.func_", GA.inputs$method), + "k", + "AIC", + "AICc", + "R2m", + "R2c", + "LL" + ) + out <- list(GA.summary = multi.GA_nG, MLPE.model = MLPE.model, + AICc.tab = AICc.tab, cd = cd.list, percent.contribution = p.cont, k = k.df) + return(out) } @@ -381,19 +404,41 @@ MS_optim <- function(CS.inputs = NULL, row.names = F, col.names = T) - file.remove(list.files(GA.inputs$Write.dir, full.names = TRUE)) + unlink(GA.inputs$Write.dir, recursive = T, force = T) k.df <- data.frame(surface = NAME, k = k) cd.list <- list(as.matrix(cd)) names(cd.list) <- NAME + AICc.tab <- data.frame(surface = NAME, + obj = multi.GA_nG@fitnessValue, + k = k, + AIC = aic, + AICc = AICc, + R2m = fit.stats[[1]], + R2c = fit.stats[[2]], + LL = LL) + + colnames(AICc.tab) <- + c( + "Surface", + paste0("obj.func_", GA.inputs$method), + "k", + "AIC", + "AICc", + "R2m", + "R2c", + "LL" + ) + out <- list(GA.summary = multi.GA_nG, MLPE.model = MLPE.model, + AICc.tab = AICc.tab, cd = cd.list, percent.contribution = p.cont, k = k.df) - + return(out) } } \ No newline at end of file diff --git a/R/MS_optim_scale.R b/R/MS_optim_scale.R index 16bcb48..ced52a2 100644 --- a/R/MS_optim_scale.R +++ b/R/MS_optim_scale.R @@ -12,12 +12,11 @@ #' @author Bill Peterman MS_optim.scale <- function(CS.inputs = NULL, gdist.inputs = NULL, - GA.inputs, - scale.surfaces = NULL) { + GA.inputs) { if (is.null(GA.inputs$scale)) { stop( - "This function should only be used if you intend to apply kernel smoothing to your resistance surfaces" + "`MS_optim.scale` should only be used if you intend to apply kernel smoothing to your resistance surfaces" ) } @@ -58,19 +57,25 @@ MS_optim.scale <- function(CS.inputs = NULL, rt <- proc.time()[3] - t1 Opt.parm <- GA.opt <- multi.GA_nG@solution + for (i in 1:GA.inputs$n.layers) { - if (GA.inputs$surface.type[i] == "cat") { - ga.p <- - GA.opt[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + 1])] - parm <- ga.p / min(ga.p) - Opt.parm[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + - 1])] <- parm - - } else { - parm <- - GA.opt[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + 1])] - Opt.parm[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + - 1])] <- parm + if (GA.inputs$surface.type[i] == "cat") { + ga.p <- + GA.opt[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + 1])] + parm <- ga.p / min(ga.p) + Opt.parm[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + + 1])] <- parm + + } else { + parm <- + GA.opt[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + 1])] + + if(length(parm == 4) & parm[4] < 0.5) { + parm[4] <- 0.000123456543210 + } + + Opt.parm[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + + 1])] <- parm } } multi.GA_nG@solution <- Opt.parm @@ -175,6 +180,7 @@ MS_optim.scale <- function(CS.inputs = NULL, n <- CS.inputs$n.Pops AICc <- (-2 * LL) + (2 * k) + ((2 * k) * (k + 1)) / (n - k - 1) + multi.GA_nG@solution[multi.GA_nG@solution == 0.000123456543210] <- 0 # Get parameter estimates MLPE.results <- MLPE.lmm_coef( @@ -203,15 +209,38 @@ MS_optim.scale <- function(CS.inputs = NULL, row.names = F, col.names = T) - file.remove(list.files(GA.inputs$Write.dir, full.names = TRUE)) + # move(list.files(GA.inputs$Write.dir, full.names = TRUE)) + unlink(GA.inputs$Write.dir, recursive = T, force = T) k.df <- data.frame(surface = NAME, k = k) cd.list <- list(as.matrix(cd)) names(cd.list) <- NAME + AICc.tab <- data.frame(surface = NAME, + obj = multi.GA_nG@fitnessValue, + k = k, + AIC = aic, + AICc = AICc, + R2m = fit.stats[[1]], + R2c = fit.stats[[2]], + LL = LL) + + colnames(AICc.tab) <- + c( + "Surface", + paste0("obj.func_", GA.inputs$method), + "k", + "AIC", + "AICc", + "R2m", + "R2c", + "LL" + ) + out <- list(GA.summary = multi.GA_nG, MLPE.model = MLPE.model, + AICc.tab = AICc.tab, cd = cd.list, percent.contribution = p.cont, k = k.df) @@ -254,12 +283,17 @@ MS_optim.scale <- function(CS.inputs = NULL, parm <- ga.p / min(ga.p) Opt.parm[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + 1])] <- parm - + } else { - parm <- - GA.opt[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + 1])] - Opt.parm[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + - 1])] <- parm + parm <- + GA.opt[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + 1])] + + if(length(parm == 4) & parm[4] < 0.5) { + parm[4] <- 0.000123456543210 + } + + Opt.parm[(GA.inputs$parm.index[i] + 1):(GA.inputs$parm.index[i + + 1])] <- parm } } multi.GA_nG@solution <- Opt.parm @@ -307,6 +341,8 @@ MS_optim.scale <- function(CS.inputs = NULL, ZZ = gdist.inputs$ZZ ) + multi.GA_nG@solution[multi.GA_nG@solution == 0.000123456543210] <- 0 + # Get parameter estimates MLPE.results <- MLPE.lmm_coef( resistance = GA.inputs$Results.dir, @@ -387,15 +423,38 @@ MS_optim.scale <- function(CS.inputs = NULL, row.names = F, col.names = T) - file.remove(list.files(GA.inputs$Write.dir, full.names = TRUE)) + # move(list.files(GA.inputs$Write.dir, full.names = TRUE)) + unlink(GA.inputs$Write.dir, recursive = T, force = T) k.df <- data.frame(surface = NAME, k = k) cd.list <- list(as.matrix(cd)) names(cd.list) <- NAME + AICc.tab <- data.frame(surface = NAME, + obj = multi.GA_nG@fitnessValue, + k = k, + AIC = aic, + AICc = AICc, + R2m = fit.stats[[1]], + R2c = fit.stats[[2]], + LL = LL) + + colnames(AICc.tab) <- + c( + "Surface", + paste0("obj.func_", GA.inputs$method), + "k", + "AIC", + "AICc", + "R2m", + "R2c", + "LL" + ) + out <- list(GA.summary = multi.GA_nG, MLPE.model = MLPE.model, + AICc.tab = AICc.tab, cd = cd.list, percent.contribution = p.cont, k = k.df) diff --git a/R/OptimFunction_Single.R b/R/OptimFunction_Single.R index 73b3b18..9f5ef45 100644 --- a/R/OptimFunction_Single.R +++ b/R/OptimFunction_Single.R @@ -217,7 +217,7 @@ Resistance.Opt_single <- if (quiet == FALSE) { cat(paste0("\t", "Iteration took ", round(rt, digits = 2), " seconds"), "\n") # cat(paste0("\t", EQ,"; ",round(SHAPE,digits=2),"; ", round(Max.SCALE,digits=2)),"\n") - cat(paste0("\t", method, " = ", round(obj.func, 4)), "\n", "\n") + cat(paste0("\t", method, " = ", round(obj.func.opt, 4)), "\n", "\n") if (!is.null(iter)) { if (GA.inputs$surface.type[iter] != "cat") { cat(paste0("\t", EQ, " | Shape = ", PARM[2], " | Max = ", PARM[3]), @@ -242,5 +242,7 @@ Resistance.Opt_single <- } } } + rm(r) + gc() return(obj.func.opt) } \ No newline at end of file diff --git a/R/OptimFunction_Single_scale.R b/R/OptimFunction_Single_scale.R index 553b858..3325071 100644 --- a/R/OptimFunction_Single_scale.R +++ b/R/OptimFunction_Single_scale.R @@ -36,6 +36,10 @@ Resistance.Opt_single.scale <- function(PARM, r <- Resistance ## Scale surface + if(PARM[4] < 0.5) { + PARM[4] <- 0.000123456543210 + } + r <- k.smooth(raster = r, sigma = PARM[4], SCALE = TRUE) @@ -48,10 +52,10 @@ Resistance.Opt_single.scale <- function(PARM, Max.SCALE <- (PARM[3]) ## If selected transformation is not in list, assign very large objective function value - + if (GA.inputs$surface.type[iter] == "cat") { stop( - "This function should only be used if you intend to apply kernel smoothing to your resistance surfaces" + "`SS_optim_scale` should only be used if you intend to apply kernel smoothing to a continuous or binary resistance surface" ) PARM <- PARM / min(PARM, na.rm = T) @@ -136,6 +140,8 @@ Resistance.Opt_single.scale <- function(PARM, # r <- reclassify(r, c(-Inf,1e-06, 1e-06,1e6,Inf,1e6)) + # CIRCUITSCAPE ------------------------------------------------------------ + if (!is.null(CS.inputs)) { writeRaster( x = r, @@ -191,7 +197,11 @@ Resistance.Opt_single.scale <- function(PARM, } ## + # gdistance --------------------------------------------------------------- + + if (!is.null(gdist.inputs)) { + cd <- Run_gdistance(gdist.inputs, r) if (method == "AIC") { @@ -236,15 +246,7 @@ Resistance.Opt_single.scale <- function(PARM, rt <- proc.time()[3] - t1 if (quiet == FALSE) { cat(paste0("\t", "Iteration took ", round(rt, digits = 2), " seconds"), "\n") - # cat(paste0("\t", EQ,"; ",round(SHAPE,digits=2),"; ", round(Max.SCALE,digits=2)),"\n") - cat(paste0("\t", method, " = ", round(obj.func, 4)), "\n", "\n") - # if (!is.null(iter)) { - # # if (GA.inputs$surface.type[iter] != "cat") { - # # cat(paste0("\t", EQ, " | Shape = ", PARM[2], " | Max = ", PARM[3]), - # # "\n", - # # "\n") - # # } - # } + cat(paste0("\t", method, " = ", round(obj.func.opt, 4)), "\n", "\n") } } else { obj.func.opt <- -99999 @@ -253,12 +255,9 @@ Resistance.Opt_single.scale <- function(PARM, if (quiet == FALSE) { cat(paste0("\t", "Iteration took ", round(rt, digits = 2), " seconds"), "\n") cat(paste0("\t", method, " = ", obj.func.opt, "\n")) - # if (!is.null(iter)) { - # if (GA.inputs$surface.type[iter] != "cat") { - # cat(paste0("EXCLUDED TRANSFORMATION", "\n", "\n")) - # } - # } } } + rm(r) + gc() return(obj.func.opt) } \ No newline at end of file diff --git a/R/Plot_trans.R b/R/Plot_trans.R index fb813df..7dec2db 100644 --- a/R/Plot_trans.R +++ b/R/Plot_trans.R @@ -1,15 +1,17 @@ #' Plot continuous surface transformation -#' +#' #' Plots a transformed continuous resistance surface against the original resistance values -#' +#' #' @param PARM Parameters to transform conintuous surface or resistance values of categorical surface. A vector of two parameters is required. The first term is the value of shape parameter (c), and the second term is the value of maximum scale parameter (b) #' @param Resistance Accepts three types of inputs. Provide either the path to the raw, untransformed resistance surface file or specify an R raster object. Alternatively, supply a vector with the minimum and manximum values (e.g., c(1,10)) #' @param transformation Transformation equation to apply. Can be provided as the name of the transformation or its numeric equivalent (see details) #' @param scale The standard deviation, in number of raster cells, to use when applying Gaussian kernel smoothing. This is the `sigma` parameter in the `kernel2dsmooth` function. (Default = NULL) #' @param print.dir Specify the directory where a .tiff of the transformation will be written (Default = NULL) +#' @param marginal.plot Logical. Should distribution plots be added to the margins of the response plot (Default = TRUE). Requires that the resistance surface be specified for `Resistance`. +#' @param marg.type Type of marginal plot to add. One of: [density, histogram, boxplot]. See \code{\link[ggExtra]{ggMarginal}} #' @param Name Name of resistance surface being transformed (optional). This will be added to the output file name. #' @return plot of transformed resistance values against original resistance values -#' @details This function will create a ggplot object and plot, so it requires \pkg{ggplot2} to be installed.\cr +#' @details This function will create a ggplot object and plot, so it requires \pkg{ggplot2} to be installed.\cr #' Equation names can be: #' \tabular{ll}{ #' \tab 1 = "Inverse-Reverse Monomolecular"\cr @@ -22,167 +24,362 @@ #' \tab 8 = "Inverse Ricker"\cr #' \tab 9 = "Distance"\cr #' } -#' +#' #' The "Distance" equation sets all cell values equal to 1. #' Because of the flexibility of the Ricker function to take a monomolecular shape (try \code{Plot.trans(PARM=c(10,100), Resistance=c(1,10), transformation="Ricker")} to see this), whenever a shape parameter >6 is selected in combination with a Ricker family transformation, the transformation reverts to a Distance transformation. In general, it seems that using a combination of intermediate Ricker and Monomolecular transformations provides the best, most flexible coverage of parameter space. This constraint has not been implemented in the \code{Plot.tans} function. -#' @usage Plot.trans(PARM, Resistance, transformation, scale, print.dir, Name) +#' @usage Plot.trans(PARM, +#' Resistance, +#' transformation, +#' scale, +#' print.dir, +#' marginal.plot, +#' marg.type, +#' Name) #' @export #' @author Bill Peterman -Plot.trans <- function(PARM,Resistance,transformation, scale = NULL, print.dir=NULL, Name="layer"){ - if(length(Resistance)==2) { +Plot.trans <- function(PARM, + Resistance, + transformation, + scale = NULL, + print.dir = NULL, + marginal.plot = TRUE, + marg.type = "histogram", + Name = "layer") { + if (length(Resistance) == 2) { + marginal.plot <- FALSE + # stop( + # "If you wish to generate marginal plots, please supply the raster surface to the `Resistance` arguement." + # ) + r <- Resistance - if(!is.null(scale)) { + if (!is.null(scale)) { zmat <- as.matrix(r) - r <- smoothie::kernel2dsmooth(zmat, - kernel.type="gauss", - nx=nrow(r), - ny=ncol(r), - sigma = scale) + + x <- spatstat::as.im(zmat) + + r <- spatstat::blur(x = x, + sigma = sigma, + normalise = TRUE, + bleed = FALSE) } - Mn=min(r) - Mx=max(r) + Mn = min(r) + Mx = max(r) NAME <- Name - } else if(class(Resistance)[1]!='RasterLayer') { - r<-raster(Resistance) + } else if (class(Resistance)[1] != 'RasterLayer') { + r <- raster(Resistance) NAME <- basename(Resistance) - NAME<-sub("^([^.]*).*", "\\1", NAME) - names(r)<-NAME - if(!is.null(scale)) { + NAME <- sub("^([^.]*).*", "\\1", NAME) + names(r) <- NAME + if (!is.null(scale)) { zmat <- as.matrix(r) - r <- smoothie::kernel2dsmooth(zmat, - kernel.type="gauss", - nx=nrow(r), - ny=ncol(r), - sigma = scale) - Mn=min(r) - Mx=max(r) + + x <- spatstat::as.im(zmat) + + r <- spatstat::blur(x = x, + sigma = sigma, + normalise = TRUE, + bleed = FALSE) + + Mn = min(r) + Mx = max(r) } else { - Mn=cellStats(r,stat='min') - Mx=cellStats(r,stat='max') + Mn = cellStats(r, stat = 'min') + Mx = cellStats(r, stat = 'max') } } else { r <- Resistance - NAME <- basename(r@file@name) - NAME<-sub("^([^.]*).*", "\\1", NAME) - names(r)<-NAME + NAME <- names(r) + # NAME <- sub("^([^.]*).*", "\\1", NAME) + # names(r) <- NAME - if(!is.null(scale)) { + if (!is.null(scale)) { zmat <- as.matrix(r) - r <- smoothie::kernel2dsmooth(zmat, - kernel.type="gauss", - nx=nrow(r), - ny=ncol(r), - sigma = scale) - Mn=min(r) - Mx=max(r) + + x <- spatstat::as.im(zmat) + + r <- spatstat::blur(x = x, + sigma = sigma, + normalise = TRUE, + bleed = FALSE) + Mn = min(r) + Mx = max(r) } else { - Mn=cellStats(r,stat='min') - Mx=cellStats(r,stat='max') + Mn = cellStats(r, stat = 'min') + Mx = cellStats(r, stat = 'max') } - } - - if(Name!="layer"){ - NAME<-Name } - # Make vector of data - original <- seq(from=Mn,to=Mx,length.out=200) - dat.t <- SCALE.vector(data=original,0,10) + if (Name != "layer") { + NAME <- Name + } - SHAPE <- PARM[[1]] - Max.SCALE <- PARM[[2]] + # Make vector data -------------------------------------------------------- - if(is.numeric(transformation)){ - equation<-get.EQ(transformation) + # * No marginal plots ----------------------------------------------------- + if(marginal.plot == FALSE) { + original <- seq(from = Mn, + to = Mx, + length.out = 200) + dat.t <- SCALE.vector(data = original, 0, 10) + + SHAPE <- PARM[[1]] + Max.SCALE <- PARM[[2]] + + if (is.numeric(transformation)) { + equation <- get.EQ(transformation) + + } else { + equation <- transformation + } + # Set equation/name combination + if (equation == "Distance") { + Trans.vec <- (dat.t * 0) + 1 + TITLE <- "Distance" + + } else if (equation == "Inverse-Reverse Monomolecular") { + SIGN <- -1 # Inverse + Trans.vec <- + SIGN * PARM[[2]] * (1 - exp(-1 * dat.t / PARM[[1]])) + SIGN # Monomolecular + Trans.vec <- + rev(SCALE.vector(Trans.vec, MIN = abs(max(Trans.vec)), MAX = abs(min(Trans.vec)))) + TITLE <- "Inverse-Reverse Monomolecular" + + } else if (equation == "Inverse Monomolecular") { + SIGN <- -1 # Inverse + Trans.vec <- + SIGN * PARM[[2]] * (1 - exp(-1 * dat.t / PARM[[1]])) + SIGN # Monomolecular + Trans.vec <- + (SCALE.vector(Trans.vec, MIN = abs(max(Trans.vec)), MAX = abs(min(Trans.vec)))) + TITLE <- "Inverse Monomolecular" + + } else if (equation == "Monomolecular") { + SIGN <- 1 + Trans.vec <- + SIGN * PARM[[2]] * (1 - exp(-1 * dat.t / PARM[[1]])) + SIGN # Monomolecular + TITLE <- "Monomolecular" + + } else if (equation == "Reverse Monomolecular") { + SIGN <- 1 + Trans.vec <- + SIGN * PARM[[2]] * (1 - exp(-1 * dat.t / PARM[[1]])) + SIGN # Monomolecular + Trans.vec <- rev(Trans.vec) # Reverse + TITLE <- "Reverse Monomolecular" + + } else if (equation == "Inverse Ricker") { + SIGN <- -1 # Inverse + Trans.vec <- + SIGN * (PARM[[2]] * dat.t * exp(-1 * dat.t / PARM[[1]])) + SIGN # Ricker + Trans.vec <- + SCALE.vector(Trans.vec, MIN = abs(max(Trans.vec)), MAX = abs(min(Trans.vec))) + TITLE <- "Inverse Ricker" + + } else if (equation == "Reverse Ricker") { + SIGN <- 1 + Trans.vec <- + SIGN * (PARM[[2]] * dat.t * exp(-1 * dat.t / PARM[[1]])) + SIGN # Ricker + Trans.vec <- rev(Trans.vec) # Reverse + TITLE <- "Reverse Ricker" + + } else if (equation == "Inverse-Reverse Ricker") { + SIGN <- -1 # Inverse + Trans.vec <- + SIGN * (PARM[[2]] * dat.t * exp(-1 * dat.t / PARM[[1]])) + SIGN # Ricker + Trans.vec <- + rev(SCALE.vector(Trans.vec, MIN = abs(max(Trans.vec)), MAX = abs(min(Trans.vec)))) # Reverse + TITLE <- "Inverse-Reverse Ricker" + + } else { + SIGN <- 1 + Trans.vec <- + SIGN * (PARM[[2]] * dat.t * exp(-1 * dat.t / PARM[[1]])) + SIGN # Ricker + TITLE <- "Ricker" + } + transformed <- Trans.vec + plot.data <- data.frame(original, transformed) + x.break <- pretty(original * 1.07) + y.break <- pretty(transformed * 1.07) + + ( + p <- ggplot(plot.data, aes(x = original, y = transformed)) + + # ggtitle(equation) + + theme_bw() + + geom_line(size = 1.5) + + xlab(expression(bold( + "Original data values" + ))) + + ylab(expression(bold( + "Transformed data values" + ))) + + theme( + # plot.title = element_text( + # lineheight = 2, + # face = "bold", + # size = 20 + # ), + legend.title = element_blank(), + legend.key = element_blank(), + axis.text.x = element_text(size = 14), + axis.text.y = element_text(size = 14), + axis.title.x = element_text(size = 16), + axis.title.y = element_text(size = 16) + ) + + scale_x_continuous(limits = c(min(original), max(original)), breaks = + x.break) + + scale_y_continuous(limits = c(min(transformed), max(transformed)), breaks = + y.break) + + removeGrid() + ) + + # * Marginal plots -------------------------------------------------------- } else { - equation<-transformation - } - # Set equation/name combination - if(equation=="Distance") { - Trans.vec <- (dat.t*0)+1 - TITLE <- "Distance" + original1 <- seq(from = Mn, + to = Mx, + length.out = 200) + dat.t1 <- SCALE.vector(data = original1, 0, 10) + type1 <- rep(2, 200) - } else if(equation=="Inverse-Reverse Monomolecular"){ - SIGN<- -1 # Inverse - Trans.vec <- SIGN*PARM[[2]]*(1-exp(-1*dat.t/PARM[[1]]))+SIGN # Monomolecular - Trans.vec <- rev(SCALE.vector(Trans.vec,MIN=abs(max(Trans.vec)),MAX=abs(min(Trans.vec)))) - TITLE <- "Inverse-Reverse Monomolecular" + original2 <- sort(as.vector(r)) + dat.t2 <- SCALE.vector(data = original2, 0, 10) + type <- rep(1, length(original2)) - } else if(equation=="Inverse Monomolecular"){ - SIGN<- -1 # Inverse - Trans.vec <- SIGN*PARM[[2]]*(1-exp(-1*dat.t/PARM[[1]]))+SIGN # Monomolecular - Trans.vec <- (SCALE.vector(Trans.vec,MIN=abs(max(Trans.vec)),MAX=abs(min(Trans.vec)))) - TITLE <- "Inverse Monomolecular" + Trans.list <- DAT <- list(dat.t2, + dat.t1 + ) - } else if(equation=="Monomolecular") { - SIGN<- 1 - Trans.vec <- SIGN*PARM[[2]]*(1-exp(-1*dat.t/PARM[[1]]))+SIGN # Monomolecular - TITLE <- "Monomolecular" - } else if(equation=="Reverse Monomolecular") { - SIGN<- 1 - Trans.vec <- SIGN*PARM[[2]]*(1-exp(-1*dat.t/PARM[[1]]))+SIGN # Monomolecular - Trans.vec <- rev(Trans.vec) # Reverse - TITLE <- "Reverse Monomolecular" + SHAPE <- PARM[[1]] + Max.SCALE <- PARM[[2]] - } else if (equation=="Inverse Ricker") { - SIGN <- -1 # Inverse - Trans.vec <- SIGN*(PARM[[2]]*dat.t*exp(-1*dat.t/PARM[[1]]))+SIGN # Ricker - Trans.vec <- SCALE.vector(Trans.vec,MIN=abs(max(Trans.vec)),MAX=abs(min(Trans.vec))) - TITLE <- "Inverse Ricker" + if (is.numeric(transformation)) { + equation <- get.EQ(transformation) + + } else { + equation <- transformation + } + # Set equation/name combination + for(i in 1:2){ + dat.t <- DAT[[i]] + if (equation == "Distance") { + Trans.vec <- (dat.t * 0) + 1 + TITLE <- "Distance" + + } else if (equation == "Inverse-Reverse Monomolecular") { + SIGN <- -1 # Inverse + Trans.vec <- + SIGN * PARM[[2]] * (1 - exp(-1 * dat.t / PARM[[1]])) + SIGN # Monomolecular + Trans.vec <- + rev(SCALE.vector(Trans.vec, MIN = abs(max(Trans.vec)), MAX = abs(min(Trans.vec)))) + TITLE <- "Inverse-Reverse Monomolecular" + + } else if (equation == "Inverse Monomolecular") { + SIGN <- -1 # Inverse + Trans.vec <- + SIGN * PARM[[2]] * (1 - exp(-1 * dat.t / PARM[[1]])) + SIGN # Monomolecular + Trans.vec <- + (SCALE.vector(Trans.vec, MIN = abs(max(Trans.vec)), MAX = abs(min(Trans.vec)))) + TITLE <- "Inverse Monomolecular" + + } else if (equation == "Monomolecular") { + SIGN <- 1 + Trans.vec <- + SIGN * PARM[[2]] * (1 - exp(-1 * dat.t / PARM[[1]])) + SIGN # Monomolecular + TITLE <- "Monomolecular" + + } else if (equation == "Reverse Monomolecular") { + SIGN <- 1 + Trans.vec <- + SIGN * PARM[[2]] * (1 - exp(-1 * dat.t / PARM[[1]])) + SIGN # Monomolecular + Trans.vec <- rev(Trans.vec) # Reverse + TITLE <- "Reverse Monomolecular" + + } else if (equation == "Inverse Ricker") { + SIGN <- -1 # Inverse + Trans.vec <- + SIGN * (PARM[[2]] * dat.t * exp(-1 * dat.t / PARM[[1]])) + SIGN # Ricker + Trans.vec <- + SCALE.vector(Trans.vec, MIN = abs(max(Trans.vec)), MAX = abs(min(Trans.vec))) + TITLE <- "Inverse Ricker" + + } else if (equation == "Reverse Ricker") { + SIGN <- 1 + Trans.vec <- + SIGN * (PARM[[2]] * dat.t * exp(-1 * dat.t / PARM[[1]])) + SIGN # Ricker + Trans.vec <- rev(Trans.vec) # Reverse + TITLE <- "Reverse Ricker" + + } else if (equation == "Inverse-Reverse Ricker") { + SIGN <- -1 # Inverse + Trans.vec <- + SIGN * (PARM[[2]] * dat.t * exp(-1 * dat.t / PARM[[1]])) + SIGN # Ricker + Trans.vec <- + rev(SCALE.vector(Trans.vec, MIN = abs(max(Trans.vec)), MAX = abs(min(Trans.vec)))) # Reverse + TITLE <- "Inverse-Reverse Ricker" + + } else { + SIGN <- 1 + Trans.vec <- + SIGN * (PARM[[2]] * dat.t * exp(-1 * dat.t / PARM[[1]])) + SIGN # Ricker + TITLE <- "Ricker" + } + Trans.list[[i]] <- Trans.vec + } # End for loop - } else if (equation=="Reverse Ricker") { - SIGN <- 1 - Trans.vec <- SIGN*(PARM[[2]]*dat.t*exp(-1*dat.t/PARM[[1]]))+SIGN # Ricker - Trans.vec <- rev(Trans.vec) # Reverse - TITLE <- "Reverse Ricker" - } else if (equation=="Inverse-Reverse Ricker") { - SIGN <- -1 # Inverse - Trans.vec <- SIGN*(PARM[[2]]*dat.t*exp(-1*dat.t/PARM[[1]]))+SIGN # Ricker - Trans.vec <- rev(SCALE.vector(Trans.vec,MIN=abs(max(Trans.vec)),MAX=abs(min(Trans.vec)))) # Reverse - TITLE <- "Inverse-Reverse Ricker" + transformed <- unlist(Trans.list) + plot.data <- data.frame(original = c(original2, original1), + transformed = transformed, + type = c(type, type1)) + x.break <- pretty(original1 * 1.07) + y.break <- pretty(transformed * 1.07) - } else { - SIGN <- 1 - Trans.vec <- SIGN*(PARM[[2]]*dat.t*exp(-1*dat.t/PARM[[1]]))+SIGN # Ricker - TITLE <- "Ricker" - } - transformed<-Trans.vec - plot.data<-data.frame(original,transformed) - x.break<-pretty(original*1.07) - y.break<-pretty(transformed*1.07) - - ( p<- ggplot(plot.data,aes(x=original,y=transformed)) + - ggtitle(equation) + + p <- ggplot(plot.data[type==1,], aes(x = original, y = transformed)) + theme_bw() + - geom_line(size=1.5) + - xlab(expression(bold("Original data values"))) + - ylab(expression(bold("Transformed data values"))) + - theme(plot.title = element_text(lineheight=2, face="bold",size=20), - legend.title = element_blank(), - legend.key = element_blank(), - axis.text.x = element_text(size=14), - axis.text.y = element_text(size=14), - axis.title.x = element_text(size=16), - axis.title.y = element_text(size=16) - # axis.line = element_line(colour = "black"), - # panel.grid.major = element_blank(), - # panel.grid.minor = element_blank(), - # panel.border = element_blank(), - # panel.background = element_blank() + geom_line(size = 1.5, color = "white") + + xlab(expression(bold( + "Original data values" + ))) + + ylab(expression(bold( + "Transformed data values" + ))) + + theme( + axis.text.x = element_text(size = 14), + axis.text.y = element_text(size = 14), + axis.title.x = element_text(size = 16), + axis.title.y = element_text(size = 16) ) + - scale_x_continuous(limits=c(min(original),max(original)),breaks=x.break) + - scale_y_continuous(limits=c(min(transformed),max(transformed)),breaks=y.break) - ) + scale_x_continuous( + breaks = x.break) + + scale_y_continuous( + breaks = y.break) + + geom_line(data = plot.data[plot.data$type==2,], + aes(x = original, y = transformed), + size = 1.5, + color = "black") + + removeGrid() + + # Add marginal plot + (p <- ggMarginal(p, + type = marg.type, + size = 8 + )) + + + } # End response plot - if(!is.null(print.dir)){ - tiff(filename=paste0(print.dir,equation,"_Transformation_",NAME ,".tif"),width=160,height=150,units="mm",res=300,compression="lzw") + if (!is.null(print.dir)) { + tiff( + filename = paste0(print.dir, equation, "_Transformation_", NAME , ".tif"), + width = 160, + height = 150, + units = "mm", + res = 300, + compression = "lzw" + ) print(p) - dev.off() + dev.off() return(p) } print(p) @@ -192,16 +389,16 @@ Plot.trans <- function(PARM,Resistance,transformation, scale = NULL, print.dir=N #### ORIGINAL #### #' #' Plot continuous surface transformation -#' #' +#' #' #' #' Plots a transformed continuous resistance surface against the original resistance values -#' #' +#' #' #' #' @param PARM Parameters to transform conintuous surface or resistance values of categorical surface. A vector of two parameters is required. The first term is the value of shape parameter (c), and the second term is the value of maximum scale parameter (b) #' #' @param Resistance Accepts three types of inputs. Provide either the path to the raw, untransformed resistance surface file or specify an R raster object. Alternatively, supply a vector with the minimum and manximum values (e.g., c(1,10)) #' #' @param transformation Transformation equation to apply. Can be provided as the name of the transformation or its numeric equivalent (see details) #' #' @param print.dir Specify the directory where a .tiff of the transformation will be written (Default = NULL) #' #' @param Name Name of resistance surface being transformed (optional). This will be added to the output file name. #' #' @return plot of transformed resistance values against original resistance values -#' #' @details This function will create a ggplot object and plot, so it requires \pkg{ggplot2} to be installed.\cr +#' #' @details This function will create a ggplot object and plot, so it requires \pkg{ggplot2} to be installed.\cr #' #' Equation names can be: #' #' \tabular{ll}{ #' #' \tab 1 = "Inverse-Reverse Monomolecular"\cr @@ -214,109 +411,109 @@ Plot.trans <- function(PARM,Resistance,transformation, scale = NULL, print.dir=N #' #' \tab 8 = "Inverse Ricker"\cr #' #' \tab 9 = "Distance"\cr #' #' } -#' #' +#' #' #' #' The "Distance" equation sets all cell values equal to 1. #' #' Because of the flexibility of the Ricker function to take a monomolecular shape (try \code{Plot.trans(PARM=c(10,100), Resistance=c(1,10), transformation="Ricker")} to see this), whenever a shape parameter >6 is selected in combination with a Ricker family transformation, the transformation reverts to a Distance transformation. In general, it seems that using a combination of intermediate Ricker and Monomolecular transformations provides the best, most flexible coverasge of parameter space. This constraint has not been implemented in the \code{Plot.tans} function. #' #' @usage Plot.trans(PARM, Resistance, transformation, print.dir, Name) #' #' @export #' #' @author Bill Peterman -#' +#' #' Plot.trans <- function(PARM,Resistance,transformation, print.dir=NULL, Name="layer"){ #' if(length(Resistance)==2) { #' r <- Resistance #' Mn=min(r) -#' Mx=max(r) +#' Mx=max(r) #' NAME <- Name #' } else if(class(Resistance)[1]!='RasterLayer') { -#' r<-raster(Resistance) +#' r<-raster(Resistance) #' NAME <- basename(Resistance) -#' NAME<-sub("^([^.]*).*", "\\1", NAME) +#' NAME<-sub("^([^.]*).*", "\\1", NAME) #' names(r)<-NAME #' Mn=cellStats(r,stat='min') -#' Mx=cellStats(r,stat='max') -#' +#' Mx=cellStats(r,stat='max') +#' #' } else { #' r <- Resistance #' Mn=cellStats(r,stat='min') #' Mx=cellStats(r,stat='max') #' NAME <- basename(r@file@name) -#' NAME<-sub("^([^.]*).*", "\\1", NAME) +#' NAME<-sub("^([^.]*).*", "\\1", NAME) #' names(r)<-NAME -#' } -#' +#' } +#' #' if(Name!="layer"){ #' NAME<-Name #' } -#' -#' # Make vector of data +#' +#' # Make vector of data #' original <- seq(from=Mn,to=Mx,length.out=200) #' dat.t <- SCALE.vector(data=original,0,10) -#' +#' #' SHAPE <- PARM[[1]] #' Max.SCALE <- PARM[[2]] -#' +#' #' if(is.numeric(transformation)){ #' equation<-get.EQ(transformation) -#' +#' #' } else { #' equation<-transformation #' } #' # Set equation/name combination #' if(equation=="Distance") { -#' Trans.vec <- (dat.t*0)+1 +#' Trans.vec <- (dat.t*0)+1 #' TITLE <- "Distance" -#' +#' #' } else if(equation=="Inverse-Reverse Monomolecular"){ #' SIGN<- -1 # Inverse #' Trans.vec <- SIGN*PARM[[2]]*(1-exp(-1*dat.t/PARM[[1]]))+SIGN # Monomolecular #' Trans.vec <- rev(SCALE.vector(Trans.vec,MIN=abs(max(Trans.vec)),MAX=abs(min(Trans.vec)))) #' TITLE <- "Inverse-Reverse Monomolecular" -#' +#' #' } else if(equation=="Inverse Monomolecular"){ #' SIGN<- -1 # Inverse #' Trans.vec <- SIGN*PARM[[2]]*(1-exp(-1*dat.t/PARM[[1]]))+SIGN # Monomolecular #' Trans.vec <- (SCALE.vector(Trans.vec,MIN=abs(max(Trans.vec)),MAX=abs(min(Trans.vec)))) #' TITLE <- "Inverse Monomolecular" -#' +#' #' } else if(equation=="Monomolecular") { #' SIGN<- 1 #' Trans.vec <- SIGN*PARM[[2]]*(1-exp(-1*dat.t/PARM[[1]]))+SIGN # Monomolecular #' TITLE <- "Monomolecular" -#' +#' #' } else if(equation=="Reverse Monomolecular") { #' SIGN<- 1 #' Trans.vec <- SIGN*PARM[[2]]*(1-exp(-1*dat.t/PARM[[1]]))+SIGN # Monomolecular #' Trans.vec <- rev(Trans.vec) # Reverse -#' TITLE <- "Reverse Monomolecular" -#' +#' TITLE <- "Reverse Monomolecular" +#' #' } else if (equation=="Inverse Ricker") { #' SIGN <- -1 # Inverse #' Trans.vec <- SIGN*(PARM[[2]]*dat.t*exp(-1*dat.t/PARM[[1]]))+SIGN # Ricker #' Trans.vec <- SCALE.vector(Trans.vec,MIN=abs(max(Trans.vec)),MAX=abs(min(Trans.vec))) -#' TITLE <- "Inverse Ricker" -#' -#' } else if (equation=="Reverse Ricker") { +#' TITLE <- "Inverse Ricker" +#' +#' } else if (equation=="Reverse Ricker") { #' SIGN <- 1 #' Trans.vec <- SIGN*(PARM[[2]]*dat.t*exp(-1*dat.t/PARM[[1]]))+SIGN # Ricker #' Trans.vec <- rev(Trans.vec) # Reverse -#' TITLE <- "Reverse Ricker" -#' -#' } else if (equation=="Inverse-Reverse Ricker") { +#' TITLE <- "Reverse Ricker" +#' +#' } else if (equation=="Inverse-Reverse Ricker") { #' SIGN <- -1 # Inverse #' Trans.vec <- SIGN*(PARM[[2]]*dat.t*exp(-1*dat.t/PARM[[1]]))+SIGN # Ricker #' Trans.vec <- rev(SCALE.vector(Trans.vec,MIN=abs(max(Trans.vec)),MAX=abs(min(Trans.vec)))) # Reverse -#' TITLE <- "Inverse-Reverse Ricker" -#' +#' TITLE <- "Inverse-Reverse Ricker" +#' #' } else { #' SIGN <- 1 #' Trans.vec <- SIGN*(PARM[[2]]*dat.t*exp(-1*dat.t/PARM[[1]]))+SIGN # Ricker #' TITLE <- "Ricker" -#' } +#' } #' transformed<-Trans.vec #' plot.data<-data.frame(original,transformed) #' x.break<-pretty(original*1.07) #' y.break<-pretty(transformed*1.07) -#' +#' #' ( p<- ggplot(plot.data,aes(x=original,y=transformed)) + #' ggtitle(equation) + #' theme_bw() + @@ -327,7 +524,7 @@ Plot.trans <- function(PARM,Resistance,transformation, scale = NULL, print.dir=N #' legend.title = element_blank(), #' legend.key = element_blank(), #' axis.text.x = element_text(size=14), -#' axis.text.y = element_text(size=14), +#' axis.text.y = element_text(size=14), #' axis.title.x = element_text(size=16), #' axis.title.y = element_text(size=16) #' # axis.line = element_line(colour = "black"), @@ -337,16 +534,16 @@ Plot.trans <- function(PARM,Resistance,transformation, scale = NULL, print.dir=N #' # panel.background = element_blank() #' ) + #' scale_x_continuous(limits=c(min(original),max(original)),breaks=x.break) + -#' scale_y_continuous(limits=c(min(transformed),max(transformed)),breaks=y.break) +#' scale_y_continuous(limits=c(min(transformed),max(transformed)),breaks=y.break) #' ) -#' +#' #' if(!is.null(print.dir)){ #' tiff(filename=paste0(print.dir,equation,"_Transformation_",NAME ,".tif"),width=160,height=150,units="mm",res=300,compression="lzw") #' print(p) -#' dev.off() +#' dev.off() #' return(p) #' } #' print(p) #' return(p) -#' +#' #' } \ No newline at end of file diff --git a/R/ResistanceGA.R b/R/ResistanceGA.R index 61dbabf..583b9d8 100644 --- a/R/ResistanceGA.R +++ b/R/ResistanceGA.R @@ -17,6 +17,8 @@ #' Official release: \url{http://www.circuitscape.org/downloads} #' #' @import raster GA lme4 ggplot2 gdistance +#' @importFrom utils combn menu +#' @importFrom ggExtra removeGrid ggMarginal #' @importFrom Matrix fac2sparse #' @importFrom plyr arrange rbind.fill ldply create_progress_bar progress_text #' @importFrom dplyr mutate group_by summarise filter tally left_join @@ -24,10 +26,10 @@ #' @importFrom MuMIn r.squaredGLMM #' @importFrom magrittr "%>%" #' @importFrom plyr arrange rbind.fill ldply -#' @importFrom smoothie kernel2dsmooth +#' @importFrom spatstat as.im blur as.matrix.im #' @importFrom grDevices dev.off tiff topo.colors #' @importFrom graphics abline filled.contour par -#' @importFrom stats AIC lm logLik qqline qqnorm resid residuals runif +#' @importFrom stats AIC lm logLik qqline qqnorm resid residuals runif as.formula sigma #' @importFrom utils file_test read.csv read.delim read.table write.table #' #' @references Please cite: diff --git a/R/Resistance_transform.R b/R/Resistance_transform.R index 536dee6..175ffdd 100644 --- a/R/Resistance_transform.R +++ b/R/Resistance_transform.R @@ -55,13 +55,12 @@ Resistance.tran <- function(transformation, if (!is.null(scale)) { zmat <- as.matrix(R) # Note that sigma is in pixels, NOT map units! - f <- smoothie::kernel2dsmooth( - zmat, - kernel.type = "gauss", - nx = nrow(R), - ny = ncol(R), - sigma = scale - ) + x <- spatstat::as.im(zmat) + + f <- spatstat::blur(x = x, + sigma = sigma, + normalise = TRUE, + bleed = FALSE) R <- R values(R) <- f } diff --git a/R/Run_gdistance.R b/R/Run_gdistance.R index 0a0583b..d50eab3 100644 --- a/R/Run_gdistance.R +++ b/R/Run_gdistance.R @@ -3,12 +3,15 @@ #' #' @param gdist.inputs Object created from running \code{\link[ResistanceGA]{CS.prep}} function #' @param r Accepts two types of inputs. Provide either the path to the resistance surface file (.asc) or specify an R RasterLayer object +#' @param scl scale the correction values (default is TRUE). No scaling should be done if the user wants to obtain absolute distance values as output. See \code{\link[gdistance]{geoCorrection}} for details #' @return A costDistance matrix object from gdistance -#' @usage Run_gdistance(gdist.inputs, r) +#' @usage Run_gdistance(gdist.inputs, r, scl) #' @export #' @author Bill Peterman -Run_gdistance <- function(gdist.inputs, r) { +Run_gdistance <- function(gdist.inputs, + r, + scl = TRUE) { if (class(r)[1] != 'RasterLayer') { r <- raster(r) } @@ -19,12 +22,12 @@ Run_gdistance <- function(gdist.inputs, r) { ) if (gdist.inputs$longlat == TRUE | gdist.inputs$directions >= 8 & gdist.inputs$method == 'costDistance') { - trC <- geoCorrection(tr, "c") + trC <- geoCorrection(tr, "c", scl = scl) ret <- costDistance(trC, gdist.inputs$samples) } if (gdist.inputs$longlat == TRUE | gdist.inputs$directions >= 8 & gdist.inputs$method == 'commuteDistance') { - trR <- geoCorrection(tr, "r") + trR <- geoCorrection(tr, "r", scl = scl) ret <- commuteDistance(trR, gdist.inputs$samples) / 1000 } diff --git a/R/SS_optim.R b/R/SS_optim.R index 2b8a7a9..b391bd0 100644 --- a/R/SS_optim.R +++ b/R/SS_optim.R @@ -14,8 +14,9 @@ #' \item A .csv file with the Maximum Likelihood Population Effects mixed effects model coefficient estimates (MLPE_coeff_Table.csv) #' \item Three summary .csv files are generated: CategoricalResults.csv, ContinuousResults.csv, & All_Results_AICc.csv. These tables contain AICc values and optimization summaries for each surface. #' } -#' All results tables are also summarized in a named list ($ContinuousResults, $CategoricalResults, $AICc, $MLPE, $MLPE.list)\cr -#' The \code{lmer} model objects stored $MLPE.list are fit using Restricted Maximum Likelihood +#' All results tables are also summarized in a named list ($ContinuousResults, $CategoricalResults, $AICc, $MLPE, $MLPE.list, $cd, $k)\cr +#' The \code{lmer} model objects stored $MLPE.list are fit using Restricted Maximum Likelihood \cr +#' $cd is a list of the optimized cost pairwise distance matrices and $k is a table of the surface names and number of parameters used to calculate AICc. These two objects can be passed to \code{\link[ResistanceGA]{Resist.boot}} to conduct a bootstrap analysis. #' @usage SS_optim(CS.inputs, gdist.inputs, GA.inputs, nlm, dist_mod, null_mod) #' @author Bill Peterman #' @export @@ -450,9 +451,6 @@ SS_optim <- function(CS.inputs = NULL, single.GA@solution[3] ) - - - colnames(RS) <- c( "Surface", @@ -685,7 +683,7 @@ SS_optim <- function(CS.inputs = NULL, names(r) <- NAME cd <- Run_gdistance(gdist.inputs, r) - save(cd, file = paste0(GA.inputs$Write.dir, NAME, ".rda")) + # save(cd, file = paste0(GA.inputs$Write.dir, NAME, ".rda")) write.table( as.matrix(cd), file = paste0(GA.inputs$Results.dir, NAME, "_", gdist.inputs$method, "_distMat.csv"), @@ -867,7 +865,7 @@ SS_optim <- function(CS.inputs = NULL, NAME <- GA.inputs$layer.names[i] cd <- Run_gdistance(gdist.inputs, r) - save(cd, file = paste0(GA.inputs$Write.dir, NAME, ".rda")) + # save(cd, file = paste0(GA.inputs$Write.dir, NAME, ".rda")) write.table( as.matrix(cd), file = paste0(GA.inputs$Results.dir, NAME, "_", gdist.inputs$method, "_distMat.csv"), @@ -954,7 +952,7 @@ SS_optim <- function(CS.inputs = NULL, NAME <- GA.inputs$layer.names[i] cd <- Run_gdistance(gdist.inputs, r) - save(cd, file = paste0(GA.inputs$Write.dir, NAME, ".rda")) + # save(cd, file = paste0(GA.inputs$Write.dir, NAME, ".rda")) write.table( as.matrix(cd), file = paste0(GA.inputs$Results.dir, NAME, "_", gdist.inputs$method, "_distMat.csv"), @@ -1414,7 +1412,7 @@ SS_optim <- function(CS.inputs = NULL, ) } - file.remove(list.files(GA.inputs$Write.dir, full.names = TRUE)) + unlink(GA.inputs$Write.dir, recursive = T, force = T) return(RESULTS) ############################################################################################################### } \ No newline at end of file diff --git a/R/SS_optim_scale.R b/R/SS_optim_scale.R index f051e71..e9aabe1 100644 --- a/R/SS_optim_scale.R +++ b/R/SS_optim_scale.R @@ -27,15 +27,21 @@ SS_optim.scale <- function(CS.inputs = NULL, null_mod = TRUE) { if (is.null(GA.inputs$scale)) { stop( - "This function should only be used if you intend to apply kernel smoothing to your resistance surfaces" + "`SS_optim.scale` should only be used if you intend to apply kernel smoothing to your resistance surfaces" + ) + } + + if (!is.null(GA.inputs$scale) & any(GA.inputs$scale.surfaces == 0)) { + stop( + "It is not currently possible to selectively scale surfaces while using 'SS_optim scale'. Either remove surfaces you do not wish to scale and/or do not specify values for `scale.surfaces` in GA.prep." ) } t1 <- proc.time()[3] - # RESULTS.cat <- list() # List to store categorical results within + RESULTS.cat <- list() # List to store categorical results within RESULTS.cont <- list() # List to store continuous results within cnt1 <- 0 - # cnt2 <- 0 + cnt2 <- 0 k.value <- GA.inputs$k.value MLPE.list <- list() cd.list <- list() @@ -62,72 +68,144 @@ SS_optim.scale <- function(CS.inputs = NULL, names(r) <- GA.inputs$layer.names[i] - single.GA <- ga( - type = "real-valued", - fitness = Resistance.Opt_single.scale, - Resistance = r, - population = GA.inputs$population, - selection = GA.inputs$selection, - pcrossover = GA.inputs$pcrossover, - pmutation = GA.inputs$pmutation, - crossover = GA.inputs$crossover, - Min.Max = GA.inputs$Min.Max, - GA.inputs = GA.inputs, - CS.inputs = CS.inputs, - min = GA.inputs$min.list[[i]], - max = GA.inputs$max.list[[i]], - popSize = GA.inputs$pop.mult * length(GA.inputs$max.list[[i]]), - maxiter = GA.inputs$maxiter, - run = GA.inputs$run, - keepBest = GA.inputs$keepBest, - elitism = GA.inputs$percent.elite, - mutation = GA.inputs$mutation, - seed = GA.inputs$seed, - iter = i, - quiet = GA.inputs$quiet - ) - - - start.vals <- single.GA@solution[-1] - EQ <- get.EQ(single.GA@solution[1]) + # Scaled optimization: CS ----------------------------------------------------- - r.tran <- - Resistance.tran( + if(GA.inputs$scale.surfaces[i] == 1) { + single.GA <- ga( + type = "real-valued", + fitness = Resistance.Opt_single.scale, + Resistance = r, + population = GA.inputs$population, + selection = GA.inputs$selection, + pcrossover = GA.inputs$pcrossover, + pmutation = GA.inputs$pmutation, + crossover = GA.inputs$crossover, + Min.Max = GA.inputs$Min.Max, + GA.inputs = GA.inputs, + CS.inputs = CS.inputs, + min = GA.inputs$min.list[[i]], + max = GA.inputs$max.list[[i]], + popSize = GA.inputs$pop.mult * length(GA.inputs$max.list[[i]]), + maxiter = GA.inputs$maxiter, + run = GA.inputs$run, + keepBest = GA.inputs$keepBest, + elitism = GA.inputs$percent.elite, + mutation = GA.inputs$mutation, + seed = GA.inputs$seed, + iter = i, + quiet = GA.inputs$quiet + ) + + start.vals <- single.GA@solution[-1] + + EQ <- get.EQ(single.GA@solution[1]) + + if(single.GA@solution[4] < 0.5) { + single.GA@solution[4] <- 0.000123456543210 + } + + r.tran <- + Resistance.tran( + transformation = single.GA@solution[1], + shape = single.GA@solution[2], + max = single.GA@solution[3], + scale = single.GA@solution[4], + r = R.orig + ) + + names(r.tran) <- GA.inputs$layer.names[i] + + Run_CS(CS.inputs, + GA.inputs, + r.tran, + EXPORT.dir = GA.inputs$Results.dir) + + Diagnostic.Plots( + resistance.mat = paste0( + GA.inputs$Results.dir, + GA.inputs$layer.names[i], + "_resistances.out" + ), + genetic.dist = CS.inputs$response, + plot.dir = GA.inputs$Plots.dir, + type = "continuous", + ID = CS.inputs$ID, + ZZ = CS.inputs$ZZ + ) + + Plot.trans( + PARM = single.GA@solution[-1], + Resistance = GA.inputs$Resistance.stack[[i]], transformation = single.GA@solution[1], - shape = single.GA@solution[2], - max = single.GA@solution[3], - scale = single.GA@solution[4], - r = R.orig + print.dir = GA.inputs$Plots.dir, + scale = single.GA@solution[4] ) - - names(r.tran) <- GA.inputs$layer.names[i] - - Run_CS(CS.inputs, - GA.inputs, - r.tran, - EXPORT.dir = GA.inputs$Results.dir) - - Diagnostic.Plots( - resistance.mat = paste0( - GA.inputs$Results.dir, - GA.inputs$layer.names[i], - "_resistances.out" - ), - genetic.dist = CS.inputs$response, - plot.dir = GA.inputs$Plots.dir, - type = "continuous", - ID = CS.inputs$ID, - ZZ = CS.inputs$ZZ - ) - - Plot.trans( - PARM = single.GA@solution[-1], - Resistance = GA.inputs$Resistance.stack[[i]], - transformation = single.GA@solution[1], - print.dir = GA.inputs$Plots.dir, - scale = single.GA@solution[4] - ) + } else { + single.GA <- ga( + type = "real-valued", + fitness = Resistance.Opt_single, + Resistance = r, + population = GA.inputs$population, + selection = GA.inputs$selection, + pcrossover = GA.inputs$pcrossover, + pmutation = GA.inputs$pmutation, + crossover = GA.inputs$crossover, + Min.Max = GA.inputs$Min.Max, + GA.inputs = GA.inputs, + CS.inputs = CS.inputs, + min = GA.inputs$min.list[[i]], + max = GA.inputs$max.list[[i]], + popSize = GA.inputs$pop.mult * length(GA.inputs$max.list[[i]]), + maxiter = GA.inputs$maxiter, + run = GA.inputs$run, + keepBest = GA.inputs$keepBest, + elitism = GA.inputs$percent.elite, + mutation = GA.inputs$mutation, + seed = GA.inputs$seed, + iter = i, + quiet = GA.inputs$quiet + ) + + start.vals <- single.GA@solution[-1] + + EQ <- get.EQ(single.GA@solution[1]) + + r.tran <- + Resistance.tran( + transformation = single.GA@solution[1], + shape = single.GA@solution[2], + max = single.GA@solution[3], + r = R.orig + ) + + names(r.tran) <- GA.inputs$layer.names[i] + + Run_CS(CS.inputs, + GA.inputs, + r.tran, + EXPORT.dir = GA.inputs$Results.dir) + + Diagnostic.Plots( + resistance.mat = paste0( + GA.inputs$Results.dir, + GA.inputs$layer.names[i], + "_resistances.out" + ), + genetic.dist = CS.inputs$response, + plot.dir = GA.inputs$Plots.dir, + type = "continuous", + ID = CS.inputs$ID, + ZZ = CS.inputs$ZZ + ) + + Plot.trans( + PARM = single.GA@solution[-1], + Resistance = GA.inputs$Resistance.stack[[i]], + transformation = single.GA@solution[1], + print.dir = GA.inputs$Plots.dir, + ) + } fit.stats <- r.squaredGLMM( MLPE.lmm( @@ -207,20 +285,42 @@ SS_optim.scale <- function(CS.inputs = NULL, AICc <- (-2 * LL) + (2 * k) + (((2 * k) * (k + 1)) / (n - k - 1)) - RS <- data.frame( - GA.inputs$layer.names[i], - single.GA@fitnessValue, - k, - aic, - AICc, - fit.stats[[1]], - fit.stats[[2]], - LL[[1]], - get.EQ(single.GA@solution[1]), - single.GA@solution[2], - single.GA@solution[3], - single.GA@solution[4] - ) + if(single.GA@solution[4] < 0.5) { + single.GA@solution[4] <- 0 + } + + if(GA.inputs$scale.surfaces[i] == 1) { + RS <- data.frame( + GA.inputs$layer.names[i], + single.GA@fitnessValue, + k, + aic, + AICc, + fit.stats[[1]], + fit.stats[[2]], + LL[[1]], + get.EQ(single.GA@solution[1]), + single.GA@solution[2], + single.GA@solution[3], + single.GA@solution[4] + ) + } else { + RS <- data.frame( + GA.inputs$layer.names[i], + single.GA@fitnessValue, + k, + aic, + AICc, + fit.stats[[1]], + fit.stats[[2]], + LL[[1]], + get.EQ(single.GA@solution[1]), + single.GA@solution[2], + single.GA@solution[3], + NA + ) + } + colnames(RS) <- c( @@ -385,182 +485,420 @@ SS_optim.scale <- function(CS.inputs = NULL, #### G-distance optimization #### if (!is.null(gdist.inputs)) { - cnt1 <- cnt1 + 1 + cnt2 <- cnt2 + 1 names(r) <- GA.inputs$layer.names[i] - single.GA <- ga( - type = "real-valued", - fitness = Resistance.Opt_single.scale, - Resistance = r, - population = GA.inputs$population, - selection = GA.inputs$selection, - pcrossover = GA.inputs$pcrossover, - pmutation = GA.inputs$pmutation, - crossover = GA.inputs$crossover, - Min.Max = GA.inputs$Min.Max, - GA.inputs = GA.inputs, - gdist.inputs = gdist.inputs, - min = GA.inputs$min.list[[i]], - max = GA.inputs$max.list[[i]], - parallel = GA.inputs$parallel, - popSize = GA.inputs$pop.mult * length(GA.inputs$max.list[[i]]), - maxiter = GA.inputs$maxiter, - run = GA.inputs$run, - keepBest = GA.inputs$keepBest, - elitism = GA.inputs$percent.elite, - mutation = GA.inputs$mutation, - seed = GA.inputs$seed, - iter = i, - quiet = GA.inputs$quiet - ) - - #!#!#! - start.vals <- single.GA@solution[-1] - - EQ <- get.EQ(single.GA@solution[1]) + # Scaled optimization: gdistance ----------------------------------------------------- - r.tran <- - Resistance.tran( - transformation = single.GA@solution[1], - shape = single.GA@solution[2], - max = single.GA@solution[3], - scale = single.GA@solution[4], - r = R.orig + if(GA.inputs$scale.surfaces[i] == 1) { + single.GA <- ga( + type = "real-valued", + fitness = Resistance.Opt_single.scale, + Resistance = r, + population = GA.inputs$population, + selection = GA.inputs$selection, + pcrossover = GA.inputs$pcrossover, + pmutation = GA.inputs$pmutation, + crossover = GA.inputs$crossover, + Min.Max = GA.inputs$Min.Max, + GA.inputs = GA.inputs, + gdist.inputs = gdist.inputs, + min = GA.inputs$min.list[[i]], + max = GA.inputs$max.list[[i]], + parallel = GA.inputs$parallel, + popSize = GA.inputs$pop.mult * length(GA.inputs$max.list[[i]]), + maxiter = GA.inputs$maxiter, + run = GA.inputs$run, + keepBest = GA.inputs$keepBest, + elitism = GA.inputs$percent.elite, + mutation = GA.inputs$mutation, + seed = GA.inputs$seed, + iter = i, + quiet = GA.inputs$quiet ) - - names(r.tran) <- GA.inputs$layer.names[i] - - NAME <- GA.inputs$layer.names[i] - - cd <- Run_gdistance(gdist.inputs, r.tran) - - save(cd, file = paste0(GA.inputs$Write.dir, NAME, ".rda")) - - write.table( - as.matrix(cd), - file = paste0(GA.inputs$Results.dir, NAME, "_", gdist.inputs$method,"_distMat.csv"), - sep = ",", - row.names = F, - col.names = F - ) - - writeRaster(r.tran, paste0(GA.inputs$Results.dir, NAME, ".asc"), overwrite = - TRUE) - - Diagnostic.Plots( - resistance.mat = cd, - genetic.dist = gdist.inputs$response, - plot.dir = GA.inputs$Plots.dir, - type = "continuous", - name = NAME, - ID = gdist.inputs$ID, - ZZ = gdist.inputs$ZZ - ) - - Plot.trans( - PARM = single.GA@solution[-1], - Resistance = GA.inputs$Resistance.stack[[i]], - transformation = single.GA@solution[1], - print.dir = GA.inputs$Plots.dir, - scale = single.GA@solution[4] - ) - - - fit.stats <- - r.squaredGLMM( - MLPE.lmm( - resistance = cd, - pairwise.genetic = gdist.inputs$response, - REML = F, + } else { # Surface not to be scaled + if (GA.inputs$surface.type[i] == 'cat') { + cnt1 <- cnt1 + 1 + names(r) <- GA.inputs$layer.names[i] + + single.GA <- ga( + type = "real-valued", + fitness = Resistance.Opt_single, + Resistance = r, + population = GA.inputs$population, + selection = GA.inputs$selection, + pcrossover = GA.inputs$pcrossover, + pmutation = GA.inputs$pmutation, + crossover = GA.inputs$crossover, + Min.Max = GA.inputs$Min.Max, + GA.inputs = GA.inputs, + gdist.inputs = gdist.inputs, + min = GA.inputs$min.list[[i]], + max = GA.inputs$max.list[[i]], + parallel = GA.inputs$parallel, + popSize = GA.inputs$pop.mult * length(GA.inputs$max.list[[i]]), + maxiter = GA.inputs$maxiter, + run = GA.inputs$run, + keepBest = GA.inputs$keepBest, + elitism = GA.inputs$percent.elite, + mutation = GA.inputs$mutation, + seed = GA.inputs$seed, + iter = i, + quiet = GA.inputs$quiet + ) + + single.GA@solution <- + single.GA@solution / min(single.GA@solution) + df <- data.frame(id = unique(r), t(single.GA@solution)) + r <- subs(r, df) + NAME <- GA.inputs$layer.names[i] + names(r) <- NAME + + cd <- Run_gdistance(gdist.inputs, r) + # save(cd, file = paste0(GA.inputs$Write.dir, NAME, ".rda")) + write.table( + as.matrix(cd), + file = paste0(GA.inputs$Results.dir, NAME, "_", gdist.inputs$method, "_distMat.csv"), + + sep = ",", + row.names = F, + col.names = F + ) + writeRaster(r, + paste0(GA.inputs$Results.dir, NAME, ".asc"), + overwrite = TRUE) + Diagnostic.Plots( + resistance.mat = cd, + genetic.dist = gdist.inputs$response, + plot.dir = GA.inputs$Plots.dir, + type = "categorical", + name = NAME, ID = gdist.inputs$ID, ZZ = gdist.inputs$ZZ ) - ) - - - aic <- - AIC( - MLPE.lmm( + + fit.stats <- r.squaredGLMM( + MLPE.lmm2( + resistance = cd, + response = gdist.inputs$response, + REML = F, + ID = gdist.inputs$ID, + ZZ = gdist.inputs$ZZ + ) + ) + + aic <- AIC( + MLPE.lmm2( + resistance = cd, + response = gdist.inputs$response, + REML = F, + ID = gdist.inputs$ID, + ZZ = gdist.inputs$ZZ + ) + ) + + LL <- logLik( + MLPE.lmm2( + resistance = cd, + response = gdist.inputs$response, + REML = F, + ID = gdist.inputs$ID, + ZZ = gdist.inputs$ZZ + ) + ) + + if (k.value == 1) { + k <- 2 + } else if (k.value == 2) { + k <- GA.inputs$parm.type$n.parm[i] + 1 + } else if (k.value == 3) { + k <- GA.inputs$parm.type$n.parm[i] + length(GA.inputs$layer.names) + 1 + } else { + k <- length(GA.inputs$layer.names[i]) + 1 + } + + k.list[[i]] <- k + names(k.list)[i] <- GA.inputs$layer.names[i] + + n <- gdist.inputs$n.Pops + AICc <- + (-2 * LL) + (2 * k) + (((2 * k) * (k + 1)) / (n - k - 1)) + + + RS <- data.frame( + GA.inputs$layer.names[i], + single.GA@fitnessValue, + k, + aic, + AICc, + fit.stats[[1]], + fit.stats[[2]], + LL[[1]], + single.GA@solution + ) + + k <- GA.inputs$parm.type$n.parm[i] + + Features <- matrix() + for (z in 1:(k)) { + feature <- paste0("Feature", z) + Features[z] <- feature + } + + colnames(RS) <- + c( + "Surface", + paste0("obj.func_", GA.inputs$method), + "k", + "AIC", + "AICc", + "R2m", + "R2c", + "LL", + Features + ) + + RESULTS.cat[[cnt1]] <- RS + + MLPE.list[[i]] <- MLPE.lmm2( resistance = cd, - pairwise.genetic = gdist.inputs$response, - REML = F, + response = gdist.inputs$response, + REML = TRUE, ID = gdist.inputs$ID, ZZ = gdist.inputs$ZZ ) - ) + + cd.list[[i]] <- as.matrix(cd) + names(cd.list)[i] <- GA.inputs$layer.names[i] + + names(MLPE.list)[i] <- GA.inputs$layer.names[i] + + } else { + # Processing of unscaled continuous surface + cnt2 <- cnt2 + 1 + r <- SCALE(r, 0, 10) + names(r) <- GA.inputs$layer.names[i] + + single.GA <- ga( + type = "real-valued", + fitness = Resistance.Opt_single, + Resistance = r, + population = GA.inputs$population, + selection = GA.inputs$selection, + pcrossover = GA.inputs$pcrossover, + pmutation = GA.inputs$pmutation, + crossover = GA.inputs$crossover, + Min.Max = GA.inputs$Min.Max, + GA.inputs = GA.inputs, + gdist.inputs = gdist.inputs, + min = GA.inputs$min.list[[i]], + max = GA.inputs$max.list[[i]], + parallel = GA.inputs$parallel, + popSize = GA.inputs$pop.mult * length(GA.inputs$max.list[[i]]), + maxiter = GA.inputs$maxiter, + run = GA.inputs$run, + keepBest = GA.inputs$keepBest, + elitism = GA.inputs$percent.elite, + mutation = GA.inputs$mutation, + seed = GA.inputs$seed, + iter = i, + quiet = GA.inputs$quiet + ) + } + } # End gdist - LL <- - logLik( - MLPE.lmm( - resistance = cd, - pairwise.genetic = gdist.inputs$response, - REML = F, - ID = gdist.inputs$ID, - ZZ = gdist.inputs$ZZ + #!#!#! + if(GA.inputs$surface.type[i] != 'cat'){ + start.vals <- single.GA@solution[-1] + + EQ <- get.EQ(single.GA@solution[1]) + + if(GA.inputs$scale.surfaces[i] == 1) { + if(single.GA@solution[4] < 0.5) { + single.GA@solution[4] <- 0.000123456543210 + } + + r.tran <- + Resistance.tran( + transformation = single.GA@solution[1], + shape = single.GA@solution[2], + max = single.GA@solution[3], + scale = single.GA@solution[4], + r = R.orig + ) + + Plot.trans( + PARM = single.GA@solution[-1], + Resistance = GA.inputs$Resistance.stack[[i]], + transformation = single.GA@solution[1], + print.dir = GA.inputs$Plots.dir, + scale = single.GA@solution[4] + ) + + } else { + r.tran <- + Resistance.tran( + transformation = single.GA@solution[1], + shape = single.GA@solution[2], + max = single.GA@solution[3], + r = R.orig + ) + + Plot.trans( + PARM = single.GA@solution[-1], + Resistance = GA.inputs$Resistance.stack[[i]], + transformation = EQ, + print.dir = GA.inputs$Plots.dir ) + } + + + names(r.tran) <- GA.inputs$layer.names[i] + + NAME <- GA.inputs$layer.names[i] + + cd <- Run_gdistance(gdist.inputs, r.tran) + + write.table( + as.matrix(cd), + file = paste0(GA.inputs$Results.dir, NAME, "_", gdist.inputs$method,"_distMat.csv"), + sep = ",", + row.names = F, + col.names = F ) - - MLPE.list[[i]] <- MLPE.lmm2( - resistance = cd, - response = gdist.inputs$response, - REML = TRUE, - ID = gdist.inputs$ID, - ZZ = gdist.inputs$ZZ - ) - - cd.list[[i]] <- as.matrix(cd) - - names(MLPE.list)[i] <- GA.inputs$layer.names[i] - names(cd.list)[i] <- GA.inputs$layer.names[i] - - if (k.value == 1) { - k <- 2 - } else if (k.value == 2) { - k <- GA.inputs$parm.type$n.parm[i] + 1 - } else if (k.value == 3) { - k <- GA.inputs$parm.type$n.parm[i] + length(GA.inputs$layer.names) + 1 - } else { - k <- length(GA.inputs$layer.names[i]) + 1 - } - - k.list[[i]] <- k - names(k.list)[i] <- GA.inputs$layer.names[i] - n <- gdist.inputs$n.Pops - AICc <- (-2 * LL) + (2 * k) + (((2 * k) * (k + 1)) / (n - k - 1)) - - RS <- data.frame( - GA.inputs$layer.names[i], - single.GA@fitnessValue, - k, - aic, - AICc, - fit.stats[[1]], - fit.stats[[2]], - LL[[1]], - get.EQ(single.GA@solution[1]), - single.GA@solution[2], - single.GA@solution[3], - single.GA@solution[4] - ) - - colnames(RS) <- - c( - "Surface", - paste0("obj.func_", GA.inputs$method), - 'k', - "AIC", - "AICc", - "R2m", - "R2c", - "LL", - "Equation", - "shape", - "max", - "scale" + + writeRaster(r.tran, paste0(GA.inputs$Results.dir, NAME, ".asc"), overwrite = + TRUE) + + Diagnostic.Plots( + resistance.mat = cd, + genetic.dist = gdist.inputs$response, + plot.dir = GA.inputs$Plots.dir, + type = "continuous", + name = NAME, + ID = gdist.inputs$ID, + ZZ = gdist.inputs$ZZ ) - RESULTS.cont[[cnt1]] <- RS - + + fit.stats <- + r.squaredGLMM( + MLPE.lmm( + resistance = cd, + pairwise.genetic = gdist.inputs$response, + REML = F, + ID = gdist.inputs$ID, + ZZ = gdist.inputs$ZZ + ) + ) + + + aic <- + AIC( + MLPE.lmm( + resistance = cd, + pairwise.genetic = gdist.inputs$response, + REML = F, + ID = gdist.inputs$ID, + ZZ = gdist.inputs$ZZ + ) + ) + + LL <- + logLik( + MLPE.lmm( + resistance = cd, + pairwise.genetic = gdist.inputs$response, + REML = F, + ID = gdist.inputs$ID, + ZZ = gdist.inputs$ZZ + ) + ) + + MLPE.list[[i]] <- MLPE.lmm2( + resistance = cd, + response = gdist.inputs$response, + REML = TRUE, + ID = gdist.inputs$ID, + ZZ = gdist.inputs$ZZ + ) + + cd.list[[i]] <- as.matrix(cd) + + names(MLPE.list)[i] <- GA.inputs$layer.names[i] + names(cd.list)[i] <- GA.inputs$layer.names[i] + + if (k.value == 1) { + k <- 2 + } else if (k.value == 2) { + k <- GA.inputs$parm.type$n.parm[i] + 1 + } else if (k.value == 3) { + k <- GA.inputs$parm.type$n.parm[i] + length(GA.inputs$layer.names) + 1 + } else { + k <- length(GA.inputs$layer.names[i]) + 1 + } + + k.list[[i]] <- k + names(k.list)[i] <- GA.inputs$layer.names[i] + n <- gdist.inputs$n.Pops + AICc <- (-2 * LL) + (2 * k) + (((2 * k) * (k + 1)) / (n - k - 1)) + + if(GA.inputs$scale.surfaces[i] == 1) { + if(single.GA@solution[4] < 0.5) { + single.GA@solution[4] <- 0 + } + + RS <- data.frame( + GA.inputs$layer.names[i], + single.GA@fitnessValue, + k, + aic, + AICc, + fit.stats[[1]], + fit.stats[[2]], + LL[[1]], + get.EQ(single.GA@solution[1]), + single.GA@solution[2], + single.GA@solution[3], + single.GA@solution[4] + ) + } else { + RS <- data.frame( + GA.inputs$layer.names[i], + single.GA@fitnessValue, + k, + aic, + AICc, + fit.stats[[1]], + fit.stats[[2]], + LL[[1]], + get.EQ(single.GA@solution[1]), + single.GA@solution[2], + single.GA@solution[3], + NA + ) + } + + + + + colnames(RS) <- + c( + "Surface", + paste0("obj.func_", GA.inputs$method), + 'k', + "AIC", + "AICc", + "R2m", + "R2c", + "LL", + "Equation", + "shape", + "max", + "scale" + ) + RESULTS.cont[[cnt2]] <- RS + } # Continuous / scaled processing if (dist_mod == TRUE) { r <- reclassify(r, c(-Inf, Inf, 1)) @@ -797,7 +1135,8 @@ SS_optim.scale <- function(CS.inputs = NULL, k = k.list ) - file.remove(list.files(GA.inputs$Write.dir, full.names = TRUE)) + # file.remove(list.files(GA.inputs$Write.dir, full.names = TRUE)) + unlink(GA.inputs$Write.dir, recursive = T, force = T) return(RESULTS) diff --git a/R/all_comb.R b/R/all_comb.R new file mode 100644 index 0000000..af1330a --- /dev/null +++ b/R/all_comb.R @@ -0,0 +1,524 @@ +#' Analyze all combinations +#' +#' A wrapper function to run both \code{\link[ResistanceGA]{SS_optim}} and \code{\link[ResistanceGA]{MS_optim}} to optimize all combinations of resistance surfaces with the Genetic Algorithms package \pkg{GA}. Following optimization, \code{\link[ResistanceGA]{Resist.boot}} is run to conduct a bootstrap analysis. This function can only be used when optimizing resistance surface with least cost path or commute distance (\pkg{gdistance}). +#' +#' @param gdist.inputs Object created from running \code{\link[ResistanceGA]{gdist.prep}} function. Defined if optimizing using gdistance +#' @param GA.inputs Object created from running \code{\link[ResistanceGA]{GA.prep}} function. Be sure that the \code{Results.dir} has been been correctly specified as "all.comb" +#' @param results.dir Directory to write and save analysis results. This should be an empty directory. Any existing files located in this directory will be deleted! +#' @param max.combination The maximum number of surfaces to include in the all combinations analysis (Default = 4). Alternatively, specify a vector with the minimum and maximum number of surfaces to combine (e.g., c(2,4). If the minimum > 1, then the single surface optimization will be skipped. +#' @param iters Number of bootstrap iterations to be conducted (Default = 1000) +#' @param sample.prop Proportion of observations to be sampled each iteration (Default = 0.75) +#' @param replicate The number of times to replicate the GA optimization process for each surface (Default = 1) +#' @param nlm Logical, if TRUE, the final step of optimization will use nlm to fine-tune parameter estimates. This may lead to overfitting in some cases. (Default = FALSE) +#' @param dist_mod Logical, if TRUE, a Distance model will be calculated and added to the output table (default = TRUE) +#' @param null_mod Logical, if TRUE, an intercept-only model will be calculated and added to the output table (Default = TRUE) + +#' @return This function optimizes resistance surfaces in isolation using \code{\link[ResistanceGA]{SS_optim}}, followed by multisurface optimization using \code{\link[ResistanceGA]{MS_optim}}, and then conducts a bootstrap analysis. +#' @usage all_comb(gdist.inputs, +#' GA.inputs, +#' results.dir, +#' max.combination = 4, +#' iters = 1000, +#' replicate = 1, +#' sample.prop = 0.75, +#' nlm = FALSE, +#' dist_mod = TRUE, +#' null_mod = TRUE) +#' @author Bill Peterman +#' @export +all_comb <- function(gdist.inputs, + GA.inputs, + results.dir, + max.combination = 4, + iters = 1000, + replicate = 1, + sample.prop = 0.75, + nlm = FALSE, + dist_mod = TRUE, + null_mod = TRUE) { + + if(!exists('results.dir')) + return(cat("ERROR: An empty results directory must be specified")) + + if(!exists('gdist.inputs')) + return(cat("ERROR: Please specify gdist.inputs")) + + if(!exists('GA.inputs')) + return(cat("ERROR: Please specify GA.inputs")) + + if(!is.null(GA.inputs$Results.dir) & + !is.null(GA.inputs$Write.dir) & + !is.null(GA.inputs$Plots.dir)) { + return(cat("ERROR: Please correctly specify the `Results.dir` as 'all.comb' when running GA.prep")) + } + + if(length(max.combination) > 2) { + return(cat("ERROR: Please specify either a single value or a vector with the minimum and maximum value")) + } + + dir.files <- list.files(results.dir) + if(length(dir.files) != 0) { + q <- yn.question(cat( + paste0("This function is about to delete all files and folders in '", results.dir, "'"), + '\n', '\n', + paste0("Do you want to proceed? Select 1 (Yes) or 2 (No), then press [Enter]"))) + + if(q == FALSE) return(cat("Function stopped")) + } + + unlink(dir(results.dir, + full.names = TRUE), + recursive = TRUE, + force = T + ) + + # Create combination list ------------------------------------------------- + mc <- max.combination + + if(length(max.combination) == 2) { + if(mc[1] == 1) { + min.combination <- 2 + max.combination <- mc[2] + ss <- TRUE + } else { + min.combination <- mc[1] + max.combination <- mc[2] + ss <- FALSE + } + } else { + min.combination <- 2 + ss <- TRUE + } + + if(max.combination > GA.inputs$n.layers) { + max.combination <- GA.inputs$n.layers + } + comb.list <- vector(mode = "list", length = (max.combination - 1)) + + + list.count <- 0 + surface.count <- 0 + for(i in min.combination:max.combination) { + list.count <- list.count + 1 + comb.list[[list.count]] <- t(combn(1:GA.inputs$n.layers, i)) + if(is.null(nrow(comb.list[[list.count]]))) { + n.comb <- 1 + } else { + n.comb <- nrow(comb.list[[list.count]]) + } + surface.count <- surface.count + n.comb + } + + all.combs <- list() + comb.names <- list() + row.index <- 0 + for(i in 1:length(comb.list)){ + combs <- comb.list[[i]] + + if(is.null(nrow(comb.list[[i]]))) { + t.combs <- 1 + } else { + t.combs <- nrow(comb.list[[i]]) + } + + for(j in 1:t.combs) { + row.index <- row.index + 1 + all.combs[[row.index]] <- combs[j,] + c.names <- GA.inputs$layer.names[combs[j,]] + comb.names[[row.index]] <- paste(c.names, collapse = ".") + } + } + + GA.input_orig <- GA.inputs + + Results <- vector(mode = 'list', length = replicate) + # Begin Replicate Loop -------------------------------------------------- + for(i in 1:replicate){ + # Skip if min combination > 1 + if(ss == FALSE) { + ss.results <- NULL + AICc.tab <- NULL + dir.create(paste0(results.dir,'rep_',i)) + } else { # Do single surface optimization + dir.create(paste0(results.dir,'rep_',i)) + dir.create(paste0(results.dir,'rep_',i, "/", "single.surface")) + dir.create(paste0(results.dir,'rep_',i, "/", "single.surface/", "Plots")) + + # Single Surface optimization ----------------------------------------------------- + if(!is.null(GA.inputs$scale)) { + + # * Single Surface: scaled -------------------------------------------------------- + + # Update GA.input directories + GA.inputs$Plots.dir <- paste0(results.dir, + 'rep_',i, + "/", + "single.surface/", + "Plots/") + + GA.inputs$Results.dir <- paste0(results.dir, + 'rep_',i, + "/", + "single.surface/") + + ss.results <- SS_optim.scale(gdist.inputs = gdist.inputs, + GA.inputs = GA.inputs, + nlm = nlm, + dist_mod = dist_mod, + null_mod = null_mod) + + AICc.tab <- ss.results$AICc + } else { + # * Single Surface -------------------------------------------------------- + + # Update GA.input directories + GA.inputs$Plots.dir <- paste0(results.dir, + 'rep_',i, + "/", + "single.surface/", + "Plots/") + + GA.inputs$Results.dir <- paste0(results.dir, + 'rep_',i, + "/", + "single.surface/") + + ss.results <- SS_optim(gdist.inputs = gdist.inputs, + GA.inputs = GA.inputs, + nlm = nlm, + dist_mod = dist_mod, + null_mod = null_mod) + + AICc.tab <- ss.results$AICc + } + } + + + # Multisurface optimization ----------------------------------------------- + + # * Multisurface: scaled ---------------------------------------------------------- + + if(!is.null(GA.inputs$scale)) { + + ms.cd <- vector(mode = 'list', + length = length(all.combs)) + + ms.k <- vector(mode = 'list', + length = length(all.combs)) + + AICc.tab_list <- vector(mode = 'list', + length = length(all.combs)) + + ms.results <- vector(mode = "list", length = length(all.combs)) + + if(is.null(ss.results)) { + n_ss.cd <- 0 + all.cd <- ms.cd + } else { + n_ss.cd <- length(ss.results$cd) + all.cd <- c(ss.results$cd, + ms.cd) + } + + + for(j in 1:length(all.combs)) { + dir.create(paste0(results.dir,'rep_',i, "/", comb.names[[j]])) + # dir.create(paste0(results.dir,'rep_',i, "/", comb.names[[j]],"/Plots")) + + # Select raster surfaces + r.vec <- 1:GA.input_orig$n.layers + drop.vec <- r.vec[!(r.vec %in% all.combs[[j]])] + asc.comb <- dropLayer(GA.input_orig$Resistance.stack, drop.vec) + + if(!is.null(GA.input_orig$inputs$select.trans)) { + s.trans <- GA.input_orig$inputs$select.trans[all.combs[[j]]] + } else { + s.trans <- GA.input_orig$inputs$select.trans + } + + if(is.null(GA.input_orig$inputs$scale)) { + ms.scale <- FALSE + } else { + ms.scale <- TRUE + } + + if(sum(GA.input_orig$inputs$scale.surfaces) == 0) { + sc.surf <- NULL + } else { + sc.surf <- GA.input_orig$inputs$scale.surfaces[all.combs[[j]]] + } + + # sc.surf <- GA.input_orig$inputs$scale.surfaces[all.combs[[j]]] + + # Update GA.input + GA.inputs <- GA.prep(ASCII.dir = asc.comb, + Results.dir = 'all.comb', + min.cat = GA.input_orig$inputs$min.cat, + max.cat = GA.input_orig$inputs$max.cat, + max.cont = GA.input_orig$inputs$max.cont, + min.scale = GA.input_orig$inputs$min.scale, + max.scale = GA.input_orig$inputs$max.scale, + cont.shape = NULL, + select.trans = s.trans, + method = GA.input_orig$inputs$method, + scale = ms.scale, + scale.surfaces = sc.surf, + k.value = GA.input_orig$inputs$k.value, + pop.mult = GA.input_orig$inputs$pop.mult, + percent.elite = GA.input_orig$inputs$percent.elite, + type = GA.input_orig$inputs$type, + pcrossover = GA.input_orig$inputs$pcrossover, + pmutation = GA.input_orig$inputs$pmutation, + maxiter = GA.input_orig$inputs$maxiter, + run = GA.input_orig$inputs$run, + keepBest = GA.input_orig$inputs$keepBest, + population = GA.input_orig$inputs$population, + selection = GA.input_orig$inputs$selection, + crossover = GA.input_orig$inputs$crossover, + mutation = GA.input_orig$inputs$mutation, + pop.size = GA.input_orig$inputs$pop.size, + parallel = GA.input_orig$inputs$parallel, + seed = GA.input_orig$inputs$seed, + quiet = GA.input_orig$inputs$quiet + ) + + # Update GA.input directories + GA.inputs$Plots.dir <- paste0(results.dir, + 'rep_',i, + "/", + comb.names[[j]], + "/") + + GA.inputs$Results.dir <- paste0(results.dir, + 'rep_',i, + "/", + comb.names[[j]], + "/") + + ms.results[[j]] <- MS_optim.scale(gdist.inputs = gdist.inputs, + GA.inputs = GA.inputs) + + all.cd[[(j + n_ss.cd)]] <- ms.results[[j]]$cd[[1]] + + ms.k[[j]] <- ms.results[[j]]$k + + AICc.tab_list[[j]] <- ms.results[[j]]$AICc.tab + } # End combination loop + + # Convert combination lists to data frames + all.k <- rbind(ss.results$k, + plyr::ldply(ms.k)) + + names(all.cd) <- all.k$surface + + if(is.null(ss.results)) { + all.AICc <- plyr::ldply(AICc.tab_list) + } else { + all.AICc <- rbind(ss.results$AICc, + plyr::ldply(AICc.tab_list)) + } + + + all.AICc <- all.AICc %>% + dplyr::arrange(., AICc) %>% + dplyr::mutate(., delta.AICc = AICc - min(AICc)) %>% + dplyr::mutate(., weight = (exp(-0.5 * delta.AICc)) / sum(exp(-0.5 * delta.AICc))) %>% + as.data.frame() + + + # * Bootstrap: scaled ----------------------------------------------------- + + obs <- gdist.inputs$n.Pops + genetic.mat <- matrix(0, obs, obs) + genetic.mat[lower.tri(genetic.mat)] <- gdist.inputs$response + + boot.results <- Resist.boot(mod.names = all.k[,1], + dist.mat = all.cd, + n.parameters = all.k[,2], + sample.prop = sample.prop, + iters = iters, + obs = obs, + genetic.mat = genetic.mat) + } else { # End scaled optimization + + # Optimization without scaling -------------------------------------------- + + # * Multisurface ---------------------------------------------------------- + + ms.cd <- vector(mode = 'list', + length = length(all.combs)) + + ms.k <- vector(mode = 'list', + length = length(all.combs)) + + AICc.tab_list <- vector(mode = 'list', + length = length(all.combs)) + + ms.results <- vector(mode = "list", length = length(all.combs)) + + n_ss.cd <- length(ss.results$cd) + all.cd <- c(ss.results$cd, + ms.cd) + + for(j in 1:length(all.combs)) { + dir.create(paste0(results.dir,'rep_',i, "/", comb.names[[j]])) + # dir.create(paste0(results.dir,'rep_',i, "/", comb.names[[j]],"/Plots")) + + # Select raster surfaces + r.vec <- 1:GA.input_orig$n.layers + drop.vec <- r.vec[!(r.vec %in% all.combs[[j]])] + asc.comb <- dropLayer(GA.input_orig$Resistance.stack, drop.vec) + + if(!is.null(GA.input_orig$inputs$select.trans)) { + s.trans <- GA.input_orig$inputs$select.trans[all.combs[[j]]] + } else { + s.trans <- GA.input_orig$inputs$select.trans + } + + if(is.null(GA.input_orig$inputs$scale)) { + ms.scale <- FALSE + } else { + ms.scale <- TRUE + } + + if(sum(GA.input_orig$inputs$scale.surfaces) == 0) { + sc.surf <- NULL + } else { + sc.surf <- GA.input_orig$inputs$scale.surfaces[all.combs[[j]]] + } + + # sc.surf <- GA.input_orig$inputs$scale.surfaces[all.combs[[j]]] + + # Update GA.input + GA.inputs <- GA.prep(ASCII.dir = asc.comb, + Results.dir = 'all.comb', + min.cat = GA.input_orig$inputs$min.cat, + max.cat = GA.input_orig$inputs$max.cat, + max.cont = GA.input_orig$inputs$max.cont, + min.scale = GA.input_orig$inputs$min.scale, + max.scale = GA.input_orig$inputs$max.scale, + cont.shape = NULL, + select.trans = s.trans, + method = GA.input_orig$inputs$method, + scale = ms.scale, + scale.surfaces = sc.surf, + k.value = GA.input_orig$inputs$k.value, + pop.mult = GA.input_orig$inputs$pop.mult, + percent.elite = GA.input_orig$inputs$percent.elite, + type = GA.input_orig$inputs$type, + pcrossover = GA.input_orig$inputs$pcrossover, + pmutation = GA.input_orig$inputs$pmutation, + maxiter = GA.input_orig$inputs$maxiter, + run = GA.input_orig$inputs$run, + keepBest = GA.input_orig$inputs$keepBest, + population = GA.input_orig$inputs$population, + selection = GA.input_orig$inputs$selection, + crossover = GA.input_orig$inputs$crossover, + mutation = GA.input_orig$inputs$mutation, + pop.size = GA.input_orig$inputs$pop.size, + parallel = GA.input_orig$inputs$parallel, + seed = GA.input_orig$inputs$seed, + quiet = GA.input_orig$inputs$quiet + ) + + # Update GA.input directories + GA.inputs$Plots.dir <- paste0(results.dir, + 'rep_',i, + "/", + comb.names[[j]], + "/") + + GA.inputs$Results.dir <- paste0(results.dir, + 'rep_',i, + "/", + comb.names[[j]], + "/") + + ms.results[[j]] <- MS_optim(gdist.inputs = gdist.inputs, + GA.inputs = GA.inputs) + + all.cd[[(j + n_ss.cd)]] <- ms.results[[j]]$cd[[1]] + + ms.k[[j]] <- ms.results[[j]]$k + + AICc.tab_list[[j]] <- ms.results[[j]]$AICc.tab + } # End combination loop + + # Convert combination lists to data frames + all.k <- rbind(ss.results$k, + plyr::ldply(ms.k)) + + names(all.cd) <- all.k$surface + + if(is.null(ss.results)) { + all.AICc <- plyr::ldply(AICc.tab_list) + } else { + all.AICc <- rbind(ss.results$AICc, + plyr::ldply(AICc.tab_list)) + } + + all.AICc <- all.AICc %>% + dplyr::arrange(., AICc) %>% + dplyr::mutate(., delta.AICc = AICc - min(AICc)) %>% + dplyr::mutate(., weight = (exp(-0.5 * delta.AICc)) / sum(exp(-0.5 * delta.AICc))) %>% + as.data.frame() + + + # * Bootstrap ----------------------------------------------------- + + obs <- gdist.inputs$n.Pops + genetic.mat <- matrix(0, obs, obs) + genetic.mat[lower.tri(genetic.mat)] <- gdist.inputs$response + + boot.results <- Resist.boot(mod.names = all.k[,1], + dist.mat = all.cd, + n.parameters = all.k[,2], + sample.prop = sample.prop, + iters = iters, + obs = obs, + genetic.mat = genetic.mat) + + } # End scaled if-else + + # Write AICc table and Boot Results to replicate results directory + write.table(x = all.AICc, + paste0(results.dir,'rep_',i,"/", + "All_Combinations_Summary.csv"), + row.names = F, + col.names = T, + sep = ",") + + write.table(x = as.data.frame(boot.results), + paste0(results.dir,'rep_',i,"/", + "Bootstrap_Results.csv"), + row.names = F, + col.names = T, + sep = ",") + + if(replicate > 1) { + Results[[i]] <- list(summary.table = all.AICc, + boot.results = boot.results, + all.k = all.k, + all.cd = all.cd, + genetic.dist_mat = genetic.mat, + ss.results = ss.results, + ms.results = ms.results + ) + names(Results)[i] <- paste0('rep_',i) + + } else { + Results <- list(summary.table = all.AICc, + boot.results = boot.results, + all.k = all.k, + all.cd = all.cd, + genetic.dist_mat = genetic.mat, + ss.results = ss.results, + ms.results = ms.results + ) + } + + } # Close replicate loop + + return(Results) + +} # End function diff --git a/R/gdistance.R b/R/gdistance.R index e9d673e..77de5bf 100644 --- a/R/gdistance.R +++ b/R/gdistance.R @@ -14,7 +14,13 @@ #' @export #' @author Bill Peterman -#' @usage gdist.prep(n.Pops, response, samples, transitionFunction, directions, longlat, method) +#' @usage gdist.prep(n.Pops, +#' response = NULL, +#' samples, +#' transitionFunction = function(x) 1 / mean(x), +#' directions = 8, +#' longlat = FALSE, +#' method = 'commuteDistance') gdist.prep <- function(n.Pops, @@ -24,8 +30,7 @@ gdist.prep <- 1 / mean(x), directions = 8, longlat = FALSE, - method = 'commuteDistance', - raster = NULL) { + method = 'commuteDistance') { if (method != 'commuteDistance') { method <- 'costDistance' diff --git a/R/internal_helper_functions.R b/R/internal_helper_functions.R index 593883a..eaa3828 100644 --- a/R/internal_helper_functions.R +++ b/R/internal_helper_functions.R @@ -230,29 +230,6 @@ read.matrix2 <- function(cs.matrix) { } -# Make to-from population list -To.From.ID <- function(POPS) { - tmp <- matrix(nrow = POPS, ncol = POPS) - dimnames(tmp) <- list(1:POPS, 1:POPS) - tmp2 <- - as.data.frame(which(row(tmp) < col(tmp), arr.ind = TRUE)) - tmp2[[2]] <- dimnames(tmp)[[2]][tmp2$col] - tmp2[[1]] <- dimnames(tmp)[[2]][tmp2$row] - colnames(tmp2) <- c("pop1", "pop2") - as.numeric(tmp2$pop1) - as.numeric(tmp2$pop2) - ID <- plyr::arrange(tmp2, as.numeric(pop1), as.numeric(pop2)) - # ID<-tmp2[with(tmp2, order(pop1, pop2)), ] - p1 <- ID[POPS - 1, 1] - p2 <- ID[POPS - 1, 2] - ID[POPS - 1, 1] <- p2 - ID[POPS - 1, 2] <- p1 - ID$pop1 <- factor(ID$pop1) - ID$pop2 <- factor(ID$pop2) - return(ID) -} - - # Create ZZ matrix for mixed effects model ZZ.mat <- function(ID) { Zl <- @@ -312,6 +289,7 @@ sv.cat <- function(levels, pop.size, min, max) { sv.cont.nG <- function(direction, pop.size, max, + min.scale, max.scale, scale = NULL) { inc <- c(1, 3) @@ -322,7 +300,7 @@ sv.cont.nG <- function(direction, if (!is.null(scale)) { cont.starts <- matrix(nrow = pop.size, ncol = 4) for (r in 1:pop.size) { - scale.parm <- runif(1, 1, max.scale) + scale.parm <- runif(1, min.scale, max.scale) if (runif(1) < .5 && direction == "Increase") { # z1<-c(sample(inc,1) z <- Increase.starts.nG(sample(inc, 1)) @@ -361,6 +339,11 @@ sv.cont.nG <- function(direction, cont.starts[r,] <- z } } + if(ncol(cont.starts) == 4) { + rs <- sample(pop.size, floor(0.25 * pop.size), replace = F) + cont.starts[rs, 4] <- 0.25 + } + cont.starts } @@ -405,8 +388,8 @@ eq.set <- function(include.list) { Please see Details of the GA.prep." ) } - } } +} get.EQ <- function(equation) { # Apply specified transformation @@ -633,4 +616,12 @@ Rev.Ricker <- function(r, parm) { rev.rast <- SCALE.vector((-1 * r), 0, 10) Ricker(rev.rast, parm) } +} + +yn.question <- function(question, add_lines_before = TRUE) { + choices <- c("Yes", "No") + if(add_lines_before) cat("------------------------\n") + the_answer <- menu(choices, title = question) + + ifelse(the_answer == 1L, TRUE, FALSE) # returns TRUE or FALSE } \ No newline at end of file diff --git a/R/k_smooth.R b/R/k_smooth.R index 0f791ec..93ce99b 100644 --- a/R/k_smooth.R +++ b/R/k_smooth.R @@ -3,7 +3,7 @@ #' Apply Gaussian kernel smoothing of specified sigma #' #' @param raster A RasterLayer object to be smoothed -#' @param sigma The standard deviation of the Gaussian smoothing parameter (see \code{\link[smoothie]{kernel2dsmooth}} documentation in the \code{smoothie} package.) +#' @param sigma The standard deviation of the Gaussian smoothing parameter (see \code{\link[spatstat]{blur}} documentation in the \code{spatstat} package.) #' @param SCALE Logical. Should the smoothed raster surface be scaled to range from 0-10 (Default = FALSE) #' @usage k.smooth (raster, sigma, SCALE) @@ -16,16 +16,29 @@ k.smooth <- function(raster, sigma, SCALE = FALSE) { - zmat <- as.matrix(raster) - f <- smoothie::kernel2dsmooth( - zmat, - kernel.type = "gauss", - nx = nrow(raster), - ny = ncol(raster), - sigma = sigma - ) + r.mat <- raster::as.matrix(raster) + x <- spatstat::as.im(r.mat) - values(raster) <- f + x.blur <- spatstat::blur(x = x, + sigma = sigma, + normalise = TRUE, + bleed = FALSE) + + + values(raster) <- spatstat::as.matrix.im(x.blur) + + + # Original using smoothie ------------------------------------------------- + # zmat <- as.matrix(raster) + # f <- smoothie::kernel2dsmooth( + # zmat, + # kernel.type = "gauss", + # nx = nrow(raster), + # ny = ncol(raster), + # sigma = sigma + # ) + # + # values(raster) <- f if(SCALE == TRUE) { diff --git a/R/lower_matrix.R b/R/lower_matrix.R index cf7e0b0..d5d1e07 100644 --- a/R/lower_matrix.R +++ b/R/lower_matrix.R @@ -1,17 +1,20 @@ #' Make a vector of the lower half of a square distance matrix -#' +#' #' This function will prepare and compile objects and commands for optimization with the GA package -#' +#' #' @param matrix Square distance matrix with no row names. #' @return A vector of the lower half of the matrix -#' +#' #' @details This is a convenience function to obtain the lower half of a matrix, which is required as input for several other functions #' @export #' @author Bill Peterman -lower<-function(matrix){ - if(is.vector(matrix)==TRUE || dim(matrix)[1]!=dim(matrix)[2]) {warning("Must provide square distance matrix with no column or row names")} - lm<-matrix[lower.tri(matrix)] +lower <- function(matrix) { + if (is.vector(matrix) == TRUE || + dim(matrix)[1] != dim(matrix)[2]) { + warning("Must provide square distance matrix with no column or row names") + } + lm <- matrix[lower.tri(matrix)] return(lm) } diff --git a/R/rga_mlpe.R b/R/rga_mlpe.R new file mode 100644 index 0000000..f5d3598 --- /dev/null +++ b/R/rga_mlpe.R @@ -0,0 +1,129 @@ +# Run Mixed effects models using formula interface +#' Run maximum likelihood population effects mixed effects model (MLPE) +#' +#' Runs MLPE as detailed by Clarke et al. (2002). This is a general function for flexibly fitting MLPE models using the standard \code{lme4} formula interface +#' +#' @param formula \code{lme4} style mixed effects model equation. +#' @param data data frame containing vectors of values from the lower half of genetic/resistance distance matrices +#' @param REML Logical. If TRUE, mixed effects model will be fit using restricted maximum likelihood. Default = FALSE +#' @param ZZ The sparse matrix object for the MLPE model. The function will automatically create this object, but it can be specified directly from the output of CS.prep or gdist.prep (Default = NULL) +#' @return A lmer object from the fitted model +#' @details An AIC value will only be returned if \code{REML = FALSE}. The random effect must be the population vector \code{pop1} generated from the function \code{\link[ResistanceGA]{To.From.ID}}. +#' @examples +#' # Create square 'distance' matrices +#' y <- matrix(rnorm(16), 4) +#' x <- matrix(rnorm(16), 4) +#' +#' # Create to-from object (4 populations sampled) +#' id <- To.From.ID(4) +#' +#' # Create data frame +#' df <- data.frame(y = lower(y), +#' x = lower(x), +#' pop = id$pop1) +#' +#' # Fit MLPE model +#' out <- mlpe_rga(formula = y ~ x + (1 | pop), +#' data = df) + +#' @export +#' @author Bill Peterman +#' @usage mlpe_rga(formula, data, REML = FALSE, ZZ = NULL) +#' @references Clarke, R. T., P. Rothery, and A. F. Raybould. 2002. Confidence limits for regression relationships between distance matrices: Estimating gene flow with distance. Journal of Agricultural, Biological, and Environmental Statistics 7:361-372. + +mlpe_rga <- + function(formula, + data, + REML = FALSE, + ZZ = NULL) { + + + if(class(formula) != 'formula') { + formula <- as.formula(formula) + } + + if(is.null(ZZ)) { + obs <- 0.5 * (sqrt((8 * nrow(data)) + 1) + 1) + ID <- To.From.ID(obs) + ZZ <- ZZ.mat(ID) + } + + # Fit model + mod <- + lFormula(formula, + data = data, + REML = REML) + mod$reTrms$Zt <- ZZ + dfun <- do.call(mkLmerDevfun, mod) + opt <- optimizeLmer(dfun) + MOD <- + (mkMerMod(environment(dfun), opt, mod$reTrms, fr = mod$fr)) + return(MOD) + } + + +# Make to-from population list +To.From.ID <- function(POPS) { + tmp <- matrix(nrow = POPS, ncol = POPS) + dimnames(tmp) <- list(1:POPS, 1:POPS) + tmp2 <- + as.data.frame(which(row(tmp) < col(tmp), arr.ind = TRUE)) + tmp2[[2]] <- dimnames(tmp)[[2]][tmp2$col] + tmp2[[1]] <- dimnames(tmp)[[2]][tmp2$row] + colnames(tmp2) <- c("pop1", "pop2") + as.numeric(tmp2$pop1) + as.numeric(tmp2$pop2) + ID <- plyr::arrange(tmp2, as.numeric(pop1), as.numeric(pop2)) + p1 <- ID[POPS - 1, 1] + p2 <- ID[POPS - 1, 2] + ID[POPS - 1, 1] <- p2 + ID[POPS - 1, 2] <- p1 + ID$pop1 <- factor(ID$pop1) + ID$pop2 <- factor(ID$pop2) + return(ID) +} + + +# Create ZZ matrix for mixed effects model +ZZ.mat <- function(ID) { + Zl <- + lapply(c("pop1", "pop2"), function(nm) + Matrix::fac2sparse(ID[[nm]], "d", drop = FALSE)) + ZZ <- Reduce("+", Zl[-1], Zl[[1]]) + return(ZZ) +} + +# Rescale function +SCALE.vector <- function(data, MIN, MAX, threshold = 1e-5) { + if (abs(MIN - MAX) < threshold) { + data[is.finite(data)] <- 0 + data + } else { + Mn = min(data) + Mx = max(data) + (MAX - MIN) / (Mx - Mn) * (data - Mx) + MAX + } +} + +# Define scaling function +# This will rescale from 1 to specified MAX +SCALE <- function(data, MIN, MAX, threshold = 1e-5) { + if (abs(MIN - MAX) < threshold) { + data[is.finite(raster::values(data))] <- 0 + data + } else { + Mn = cellStats(data, stat = 'min') + Mx = cellStats(data, stat = 'max') + (MAX - MIN) / (Mx - Mn) * (data - Mx) + MAX + } +} + + +lower <- function(matrix) { + if (is.vector(matrix) == TRUE || + dim(matrix)[1] != dim(matrix)[2]) { + warning("Must provide square distance matrix with no column or row names") + } + lm <- matrix[lower.tri(matrix)] + return(lm) +} \ No newline at end of file diff --git a/R/to_from.R b/R/to_from.R new file mode 100644 index 0000000..da3797a --- /dev/null +++ b/R/to_from.R @@ -0,0 +1,32 @@ +# Make to-from population list +#' Convenience function to make to-from object needed for specifying MLPE random effects +#' +#' @param POPS Provide an integer value representing the number of individuals or populations that were sampled +#' @return A data frame with two columns +#' @details This object indicates the population pairs that distance values were calculated between. Note: Distance values must be taken from the lower half of a distance matrix. You can use \code{\link[ResistanceGA]{lower}} to obtain these values. When specifying a MLPE model using \code{\link[ResistanceGA]{mlpe_rga}}, \code{pop1} from \code{To.From.ID} must be added to the data frame containing model data, and must be specified as the random effect. +#' +#' @examples To.From.ID(3) + +#' @export +#' @author Bill Peterman +#' @usage To.From.ID(POPS) + +To.From.ID <- function(POPS) { + tmp <- matrix(nrow = POPS, ncol = POPS) + dimnames(tmp) <- list(1:POPS, 1:POPS) + tmp2 <- + as.data.frame(which(row(tmp) < col(tmp), arr.ind = TRUE)) + tmp2[[2]] <- dimnames(tmp)[[2]][tmp2$col] + tmp2[[1]] <- dimnames(tmp)[[2]][tmp2$row] + colnames(tmp2) <- c("pop1", "pop2") + as.numeric(tmp2$pop1) + as.numeric(tmp2$pop2) + ID <- plyr::arrange(tmp2, as.numeric(pop1), as.numeric(pop2)) + p1 <- ID[POPS - 1, 1] + p2 <- ID[POPS - 1, 2] + ID[POPS - 1, 1] <- p2 + ID[POPS - 1, 2] <- p1 + ID$pop1 <- factor(ID$pop1) + ID$pop2 <- factor(ID$pop2) + return(ID) +} \ No newline at end of file diff --git a/man/Combine_Surfaces.Rd b/man/Combine_Surfaces.Rd index e46906a..c523ff1 100644 --- a/man/Combine_Surfaces.Rd +++ b/man/Combine_Surfaces.Rd @@ -4,7 +4,14 @@ \alias{Combine_Surfaces} \title{Combine multiple resistance surfaces together} \usage{ -Combine_Surfaces(PARM, CS.inputs, gdist.inputs, GA.inputs, out, File.name, rescale, p.contribution) +Combine_Surfaces(PARM, + CS.inputs, + gdist.inputs, + GA.inputs, + out, + File.name, + rescale, + p.contribution) } \arguments{ \item{PARM}{Parameters to transform conintuous surface or resistance values of categorical surface. Requires a vector with parameters specified in the order of resistance surfaces} diff --git a/man/Combine_Surfaces.scale.Rd b/man/Combine_Surfaces.scale.Rd index 44f874f..792e785 100644 --- a/man/Combine_Surfaces.scale.Rd +++ b/man/Combine_Surfaces.scale.Rd @@ -4,7 +4,14 @@ \alias{Combine_Surfaces.scale} \title{Combine multiple resistance surfaces together} \usage{ -Combine_Surfaces.scale(PARM, CS.inputs, gdist.inputs, GA.inputs, out, File.name, rescale, p.contribution) +Combine_Surfaces.scale(PARM, + CS.inputs, + gdist.inputs, + GA.inputs, + out, + File.name, + rescale, + p.contribution) } \arguments{ \item{PARM}{Parameters to transform conintuous surface or resistance values of categorical surface. Requires a vector with parameters specified in the order of resistance surfaces} diff --git a/man/GA.prep.Rd b/man/GA.prep.Rd index a00a109..445889c 100644 --- a/man/GA.prep.Rd +++ b/man/GA.prep.Rd @@ -37,7 +37,7 @@ quiet = FALSE) \arguments{ \item{ASCII.dir}{Directory containing all raster objects to optimized. If optimizing using least cost paths, a RasterStack or RasterLayer object can be specified.} -\item{Results.dir}{If a RasterStack is provided in place of a directory containing .asc files for ASCII.dir, then a directory to export optimization results must be specified. It is critical that there are NO SPACES in the directory, as this will cause the function to fail.} +\item{Results.dir}{If a RasterStack is provided in place of a directory containing .asc files for ASCII.dir, then a directory to export optimization results must be specified. It is critical that there are NO SPACES in the directory, as this will cause the function to fail. If using the \code{\link[ResistanceGA]{all_comb}} function, specify \code{Results.dir} as "all.comb".} \item{min.cat}{The minimum value to be assessed during optimization of of categorical resistance surfaces (Default = 1e-04)} @@ -45,11 +45,11 @@ quiet = FALSE) \item{max.cont}{The maximum value to be assessed during optimization of of continuous resistance surfaces (Default = 2500)} -\item{min.scale}{The minimum scaling parameter value to be assessed during optimization of resistance surfaces with kernel smoothing (Default = 1)} +\item{min.scale}{The minimum scaling parameter value to be assessed during optimization of resistance surfaces with kernel smoothing (Default = 0.01). See details} -\item{max.scale}{The maximum scaling parameter value to be assessed during optimization of resistance surfaces with kernel smoothing (Default = 0.25 * nrows in the raster surface)} +\item{max.scale}{The maximum scaling parameter value to be assessed during optimization of resistance surfaces with kernel smoothing (Default = 0.1 * maximum dimension of the raster surface)} -\item{cont.shape}{A vector of hypothesized relationships that each continuous resistance surface will have in relation to the genetic distance reposnse (Default = NULL; see details)} +\item{cont.shape}{A vector of hypothesized relationships that each continuous resistance surface will have in relation to the genetic distance response (Default = NULL; see details)} \item{select.trans}{Option to specify which transformations are applied to continuous surfaces. Must be provided as a list. "A" = All, "M" = Monomolecular only, "R" = Ricker only. See Details.} @@ -99,7 +99,7 @@ quiet = FALSE) \item{seed}{Integer random number seed to replicate \code{ga} optimization} -\item{quiet}{Logical. If TRUE, AICc and step run time will not be printed to the screen after each step. Only \code{ga} summary information will be printed following each iteration. Default = FALSE} +\item{quiet}{Logical. If TRUE, the objective function and step run time will not be printed to the screen after each step. Only \code{ga} summary information will be printed following each iteration. (Default = FALSE)} } \value{ An R object that is a required input into optimization functions @@ -114,6 +114,8 @@ When \code{scale = TRUE}, the standard deviation of the Gaussian kernel smoothin \code{scale.surfaces} can be used to specify which surfaces to apply kernel smoothing to during multisurface optimization. For example, \code{scale.surfaces = c(1, 0, 1)} will result in the first and third surfaces being optimized with a kernel smoothing function, while the second surface will not be scaled. The order of surfaces will match either the order of the raster stack, or alphabetical order when reading in from a directory. +\code{min.scale} defaults to a minimum of 0.01. During optimization, whenever the scaling factor (sigma) is less than 0.5, ResistanceGA will not apply scaling. In this way, it is possible for a surface to not be scaled. + The Default for \code{k.value} is 2, which sets k equal to the number of parameters optimized, plus 1 for the intercept term. Prior to version 3.0-0, \code{k.value} could not be specified by the user and followed setting 2, such that k was equal to the number of parameters optimized plus the intercept term. \code{cont.shape} can take values of "Increase", "Decrease", or "Peaked". If you believe a resistance surface is related to your response in a particular way, specifying this here may decrease the time to optimization. \code{cont.shape} is used to generate an initial set of parameter values to test during optimization. If specified, a greater proportion of the starting values will include your believed relatiosnship. If unspecified (the Default), a completely random set of starting values will be generated. diff --git a/man/MLPE.lmm.Rd b/man/MLPE.lmm.Rd index 79f5ffc..e39d9b3 100644 --- a/man/MLPE.lmm.Rd +++ b/man/MLPE.lmm.Rd @@ -4,7 +4,12 @@ \alias{MLPE.lmm} \title{Run maximum likelihood population effects mixed effects model (MLPE)} \usage{ -MLPE.lmm(resistance, pairwise.genetic, REML, ID, ZZ) +MLPE.lmm(resistance, + pairwise.genetic, + REML = FALSE, + ID = NULL, + ZZ = NULL, + scale = TRUE) } \arguments{ \item{resistance}{Path to pairwise resistance distance matrix (resistances.out) from CS results. Alternatively, provide the pairwise resistances created from optimizing with `gdistance` (result of Run_gdistance).} @@ -23,7 +28,7 @@ MLPE.lmm(resistance, pairwise.genetic, REML, ID, ZZ) A lmer object from the fitted model } \description{ -Runs MLPE as detailed by Clarke et al. (2002). This function will run the model and return glmer object +Runs MLPE as detailed by Clarke et al. (2002). This function will run the model and return lmer object } \details{ An AIC value will only be returned if \code{REML = FALSE} diff --git a/man/Plot.trans.Rd b/man/Plot.trans.Rd index b5b18d7..cc2fb72 100644 --- a/man/Plot.trans.Rd +++ b/man/Plot.trans.Rd @@ -4,7 +4,14 @@ \alias{Plot.trans} \title{Plot continuous surface transformation} \usage{ -Plot.trans(PARM, Resistance, transformation, scale, print.dir, Name) +Plot.trans(PARM, + Resistance, + transformation, + scale, + print.dir, + marginal.plot, + marg.type, + Name) } \arguments{ \item{PARM}{Parameters to transform conintuous surface or resistance values of categorical surface. A vector of two parameters is required. The first term is the value of shape parameter (c), and the second term is the value of maximum scale parameter (b)} @@ -17,6 +24,10 @@ Plot.trans(PARM, Resistance, transformation, scale, print.dir, Name) \item{print.dir}{Specify the directory where a .tiff of the transformation will be written (Default = NULL)} +\item{marginal.plot}{Logical. Should distribution plots be added to the margins of the response plot (Default = TRUE). Requires that the resistance surface be specified for `Resistance`.} + +\item{marg.type}{Type of marginal plot to add. One of: [density, histogram, boxplot]. See \code{\link[ggExtra]{ggMarginal}}} + \item{Name}{Name of resistance surface being transformed (optional). This will be added to the output file name.} } \value{ @@ -26,7 +37,7 @@ plot of transformed resistance values against original resistance values Plots a transformed continuous resistance surface against the original resistance values } \details{ -This function will create a ggplot object and plot, so it requires \pkg{ggplot2} to be installed.\cr +This function will create a ggplot object and plot, so it requires \pkg{ggplot2} to be installed.\cr Equation names can be: \tabular{ll}{ \tab 1 = "Inverse-Reverse Monomolecular"\cr @@ -39,7 +50,7 @@ Equation names can be: \tab 8 = "Inverse Ricker"\cr \tab 9 = "Distance"\cr } - + The "Distance" equation sets all cell values equal to 1. Because of the flexibility of the Ricker function to take a monomolecular shape (try \code{Plot.trans(PARM=c(10,100), Resistance=c(1,10), transformation="Ricker")} to see this), whenever a shape parameter >6 is selected in combination with a Ricker family transformation, the transformation reverts to a Distance transformation. In general, it seems that using a combination of intermediate Ricker and Monomolecular transformations provides the best, most flexible coverage of parameter space. This constraint has not been implemented in the \code{Plot.tans} function. } diff --git a/man/Resist.boot.Rd b/man/Resist.boot.Rd index 02ed7f4..b97e599 100644 --- a/man/Resist.boot.Rd +++ b/man/Resist.boot.Rd @@ -4,7 +4,13 @@ \alias{Resist.boot} \title{Run bootstrap on optimized resistance surfaces} \usage{ -Resist.boot(mod.names, dist.mat, n.parameters, sample.prop, iters, obs, genetic.mat) +Resist.boot(mod.names, + dist.mat, + n.parameters, + sample.prop, + iters, obs, + rank.method, + genetic.mat) } \arguments{ \item{mod.names}{A vector of the model names to be assessed} @@ -19,16 +25,18 @@ Resist.boot(mod.names, dist.mat, n.parameters, sample.prop, iters, obs, genetic. \item{obs}{Total number of observations (populations or individuals) in your original analysis} +\item{rank.method}{What metric should be used to rank models during bootstrap analysis? c('AIC', 'AICc', 'R2', 'LL'). Default = 'AICc'} + \item{genetic.mat}{Genetic distance matrix without row or column names.} } \value{ -A data frame reporting the average model weight, average rank, number of time a model was the top model in teh set, and the frequency a model was best. +A data frame reporting the average model weight, average rank, number of time a model was the top model in the set, and the frequency a model was best. } \description{ This is a 'pseudo' bootstrap procedure that will subsample a specified proportion of the sample locations/individuals, and will refit the MLPE model using the previously optimized resistance surface. } \details{ -This is a 'pseudo'bootstrap procedure that subsamples distance and genetic matrices, refits the MLPE model for each surface, and assesses the model AICc. AICc is calculated based on the number of parameters specified. +This is a 'pseudo'bootstrap procedure that subsamples distance and genetic matrices, refits the MLPE model for each surface. AICc is calculated based on the number of parameters specified. Ranking of models during the bootstrap analysis is based on the specified \code{rank.method}, which defaults to 'AICc'. The objective of this procedure is to identify the surfaces that is top ranked across all bootstrap iterations. } \author{ Bill Peterman diff --git a/man/Run_gdistance.Rd b/man/Run_gdistance.Rd index 5ade6db..f74c7cf 100644 --- a/man/Run_gdistance.Rd +++ b/man/Run_gdistance.Rd @@ -5,12 +5,14 @@ \title{Get cost distance using gdistance Execute gdistance} \usage{ -Run_gdistance(gdist.inputs, r) +Run_gdistance(gdist.inputs, r, scl) } \arguments{ \item{gdist.inputs}{Object created from running \code{\link[ResistanceGA]{CS.prep}} function} \item{r}{Accepts two types of inputs. Provide either the path to the resistance surface file (.asc) or specify an R RasterLayer object} + +\item{scl}{scale the correction values (default is TRUE). No scaling should be done if the user wants to obtain absolute distance values as output. See \code{\link[gdistance]{geoCorrection}} for details} } \value{ A costDistance matrix object from gdistance diff --git a/man/SS_optim.Rd b/man/SS_optim.Rd index c7c6088..8cdfc95 100644 --- a/man/SS_optim.Rd +++ b/man/SS_optim.Rd @@ -26,8 +26,9 @@ This function optimizes resistance surfaces in isolation. Following optimization \item A .csv file with the Maximum Likelihood Population Effects mixed effects model coefficient estimates (MLPE_coeff_Table.csv) \item Three summary .csv files are generated: CategoricalResults.csv, ContinuousResults.csv, & All_Results_AICc.csv. These tables contain AICc values and optimization summaries for each surface. } -All results tables are also summarized in a named list ($ContinuousResults, $CategoricalResults, $AICc, $MLPE, $MLPE.list)\cr -The \code{lmer} model objects stored $MLPE.list are fit using Restricted Maximum Likelihood +All results tables are also summarized in a named list ($ContinuousResults, $CategoricalResults, $AICc, $MLPE, $MLPE.list, $cd, $k)\cr +The \code{lmer} model objects stored $MLPE.list are fit using Restricted Maximum Likelihood \cr +$cd is a list of the optimized cost pairwise distance matrices and $k is a table of the surface names and number of parameters used to calculate AICc. These two objects can be passed to \code{\link[ResistanceGA]{Resist.boot}} to conduct a bootstrap analysis. } \description{ Optimize all surfaces contained in a directory using a genetic algorithm executed with the \code{\link[GA]{ga}} function in the Genetic Algorithms package \pkg{GA} diff --git a/man/To.From.ID.Rd b/man/To.From.ID.Rd new file mode 100644 index 0000000..0d14b3a --- /dev/null +++ b/man/To.From.ID.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/to_from.R +\name{To.From.ID} +\alias{To.From.ID} +\title{Convenience function to make to-from object needed for specifying MLPE random effects} +\usage{ +To.From.ID(POPS) +} +\arguments{ +\item{POPS}{Provide an integer value representing the number of individuals or populations that were sampled} +} +\value{ +A data frame with two columns +} +\description{ +Convenience function to make to-from object needed for specifying MLPE random effects +} +\details{ +This object indicates the population pairs that distance values were calculated between. Note: Distance values must be taken from the lower half of a distance matrix. You can use \code{\link[ResistanceGA]{lower}} to obtain these values. When specifying a MLPE model using \code{\link[ResistanceGA]{mlpe_rga}}, \code{pop1} from \code{To.From.ID} must be added to the data frame containing model data, and must be specified as the random effect. +} +\examples{ +To.From.ID(3) +} +\author{ +Bill Peterman +} diff --git a/man/all_comb.Rd b/man/all_comb.Rd new file mode 100644 index 0000000..2ea0003 --- /dev/null +++ b/man/all_comb.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/all_comb.R +\name{all_comb} +\alias{all_comb} +\title{Analyze all combinations} +\usage{ +all_comb(gdist.inputs, + GA.inputs, + results.dir, + max.combination = 4, + iters = 1000, + replicate = 1, + sample.prop = 0.75, + nlm = FALSE, + dist_mod = TRUE, + null_mod = TRUE) +} +\arguments{ +\item{gdist.inputs}{Object created from running \code{\link[ResistanceGA]{gdist.prep}} function. Defined if optimizing using gdistance} + +\item{GA.inputs}{Object created from running \code{\link[ResistanceGA]{GA.prep}} function. Be sure that the \code{Results.dir} has been been correctly specified as "all.comb"} + +\item{results.dir}{Directory to write and save analysis results. This should be an empty directory. Any existing files located in this directory will be deleted!} + +\item{max.combination}{The maximum number of surfaces to include in the all combinations analysis (Default = 4). Alternatively, specify a vector with the minimum and maximum number of surfaces to combine (e.g., c(2,4). If the minimum > 1, then the single surface optimization will be skipped.} + +\item{iters}{Number of bootstrap iterations to be conducted (Default = 1000)} + +\item{replicate}{The number of times to replicate the GA optimization process for each surface (Default = 1)} + +\item{sample.prop}{Proportion of observations to be sampled each iteration (Default = 0.75)} + +\item{nlm}{Logical, if TRUE, the final step of optimization will use nlm to fine-tune parameter estimates. This may lead to overfitting in some cases. (Default = FALSE)} + +\item{dist_mod}{Logical, if TRUE, a Distance model will be calculated and added to the output table (default = TRUE)} + +\item{null_mod}{Logical, if TRUE, an intercept-only model will be calculated and added to the output table (Default = TRUE)} +} +\value{ +This function optimizes resistance surfaces in isolation using \code{\link[ResistanceGA]{SS_optim}}, followed by multisurface optimization using \code{\link[ResistanceGA]{MS_optim}}, and then conducts a bootstrap analysis. +} +\description{ +A wrapper function to run both \code{\link[ResistanceGA]{SS_optim}} and \code{\link[ResistanceGA]{MS_optim}} to optimize all combinations of resistance surfaces with the Genetic Algorithms package \pkg{GA}. Following optimization, \code{\link[ResistanceGA]{Resist.boot}} is run to conduct a bootstrap analysis. This function can only be used when optimizing resistance surface with least cost path or commute distance (\pkg{gdistance}). +} +\author{ +Bill Peterman +} diff --git a/man/gdist.prep.Rd b/man/gdist.prep.Rd index 9f980f7..34a1084 100644 --- a/man/gdist.prep.Rd +++ b/man/gdist.prep.Rd @@ -4,7 +4,13 @@ \alias{gdist.prep} \title{Prepare data for optimization using \code{gdistance}} \usage{ -gdist.prep(n.Pops, response, samples, transitionFunction, directions, longlat, method) +gdist.prep(n.Pops, + response = NULL, + samples, + transitionFunction = function(x) 1 / mean(x), + directions = 8, + longlat = FALSE, + method = 'commuteDistance') } \arguments{ \item{n.Pops}{The number of populations that are being assessed} diff --git a/man/k.smooth.Rd b/man/k.smooth.Rd index bfded9d..50c2761 100644 --- a/man/k.smooth.Rd +++ b/man/k.smooth.Rd @@ -9,7 +9,7 @@ k.smooth (raster, sigma, SCALE) \arguments{ \item{raster}{A RasterLayer object to be smoothed} -\item{sigma}{The standard deviation of the Gaussian smoothing parameter (see \code{\link[smoothie]{kernel2dsmooth}} documentation in the \code{smoothie} package.)} +\item{sigma}{The standard deviation of the Gaussian smoothing parameter (see \code{\link[spatstat]{blur}} documentation in the \code{spatstat} package.)} \item{SCALE}{Logical. Should the smoothed raster surface be scaled to range from 0-10 (Default = FALSE)} } diff --git a/man/mlpe_rga.Rd b/man/mlpe_rga.Rd new file mode 100644 index 0000000..ab7593d --- /dev/null +++ b/man/mlpe_rga.Rd @@ -0,0 +1,50 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/rga_mlpe.R +\name{mlpe_rga} +\alias{mlpe_rga} +\title{Run maximum likelihood population effects mixed effects model (MLPE)} +\usage{ +mlpe_rga(formula, data, REML = FALSE, ZZ = NULL) +} +\arguments{ +\item{formula}{\code{lme4} style mixed effects model equation.} + +\item{data}{data frame containing vectors of values from the lower half of genetic/resistance distance matrices} + +\item{REML}{Logical. If TRUE, mixed effects model will be fit using restricted maximum likelihood. Default = FALSE} + +\item{ZZ}{The sparse matrix object for the MLPE model. The function will automatically create this object, but it can be specified directly from the output of CS.prep or gdist.prep (Default = NULL)} +} +\value{ +A lmer object from the fitted model +} +\description{ +Runs MLPE as detailed by Clarke et al. (2002). This is a general function for flexibly fitting MLPE models using the standard \code{lme4} formula interface +} +\details{ +An AIC value will only be returned if \code{REML = FALSE}. The random effect must be the population vector \code{pop1} generated from the function \code{\link[ResistanceGA]{To.From.ID}}. +} +\examples{ + +# Create square 'distance' matrices +y <- matrix(rnorm(16), 4) +x <- matrix(rnorm(16), 4) + +# Create to-from object (4 populations sampled) +id <- To.From.ID(4) + +# Create data frame +df <- data.frame(y = lower(y), + x = lower(x), + pop = id$pop1) + +# Fit MLPE model +out <- mlpe_rga(formula = y ~ x + (1 | pop), + data = df) +} +\references{ +Clarke, R. T., P. Rothery, and A. F. Raybould. 2002. Confidence limits for regression relationships between distance matrices: Estimating gene flow with distance. Journal of Agricultural, Biological, and Environmental Statistics 7:361-372. +} +\author{ +Bill Peterman +} diff --git a/vignettes/ResistanceGA.R b/vignettes/ResistanceGA.R index 4eeba01..cc6f753 100644 --- a/vignettes/ResistanceGA.R +++ b/vignettes/ResistanceGA.R @@ -1,8 +1,8 @@ ## ----setup, include=FALSE------------------------------------------------ - knitr::opts_chunk$set(echo = TRUE) +knitr::opts_chunk$set(echo = TRUE) ## ---- echo = FALSE------------------------------------------------------- - library(ggplot2, warn.conflicts = F, quietly = T) +library(ggplot2, warn.conflicts = F, quietly = T) library(raster, warn.conflicts = F, quietly = T) ## ----install.package, eval=FALSE----------------------------------------- @@ -377,17 +377,9 @@ plot(sample.locales, pch=16, col="blue", add=TRUE) # # # The 'true' resistance surface will be the composite surface # # Combine resistance surfaces, omitting the feature surface -<<<<<<< HEAD # # Use an Inverse-Reverse Monomolecular transformation of the continuous surface # # Inverse-Reverse Monomolecular = 1 # PARM <- c(1, 250, 100, 1, 1.5, 150) -======= -# # Use an Reverse Monomolecular transformation of the continuous surface -# # Reverse Monomolecular = 5 -# PARM <- c(1, 250, 100, 1, 1.5, 150) -# # PARM <- c(1, 250, 100, 1, 1.5, 150) # GOOD ->>>>>>> 18eb1de80bf9b72ab953f640889493ec0c87b31e -# # # # Setting `p.contribution = TRUE` to see how each surface # # contributes to the total resistance of the composite surface diff --git a/vignettes/ResistanceGA.Rmd b/vignettes/ResistanceGA.Rmd index 76a164c..bb3b8f3 100644 --- a/vignettes/ResistanceGA.Rmd +++ b/vignettes/ResistanceGA.Rmd @@ -7,15 +7,15 @@ output: html_document: rmarkdown::html_vignette vignette: | %\VignetteIndexEntry{ResistanceGA User Guide} - %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown} --- ```{r setup, include=FALSE} - knitr::opts_chunk$set(echo = TRUE) +knitr::opts_chunk$set(echo = TRUE) ``` - -```{r, echo = FALSE} - library(ggplot2, warn.conflicts = F, quietly = T) + +```{r, echo = FALSE, warning=F} +library(ggplot2, warn.conflicts = F, quietly = T) library(raster, warn.conflicts = F, quietly = T) ``` @@ -67,7 +67,7 @@ Setup ### Install necessary software and packages If you want to optimize using CIRCUITSCAPE, this package requires that you have [CIRCUITSCAPE v4.0](http://www.circuitscape.org/downloads "CS Downloads") or higher installed on your Windows machine. At this point in time, `ResistanceGA` can only execute CIRCUITSCAPE on *WINDOWS* operating systems. You will also need to have [R >= v3.0](http://www.r-project.org/ "R downloads") installed. I would highly recommend installing [R studio](https://www.rstudio.com/ide/download/ "R Studio download") when working with R. -__Because it has now been determined that `commuteDistance` is equivalent to CIRCUITSCAPE, it is highly recommended to to conduct all optimizations using this approach, run in parallel. Doing so allows for all processes to be conducted within R, and will not require execution of software outside of the R environment. If final current flow maps are desired, these can still be generated using the optimized resistance surface using CIRCUITSCAPE.__ +__Because it has now been determined that `commuteDistance` is equivalent to CIRCUITSCAPE, it is highly recommended to conduct all optimizations using this approach, run in parallel. Doing so allows for all processes to be conducted within R, and will not require execution of software outside of the R environment. If final current flow maps are desired, these can still be generated using the optimized resistance surface using CIRCUITSCAPE.__ ### Running on Linux @@ -184,7 +184,7 @@ GA.inputs <- GA.prep(ASCII.dir = write.dir, max.cat = 500, max.cont = 500, select.trans = "M", - method = "LL" + method = "LL", seed = 555) CS.inputs <- CS.prep(n.Pops = length(sample.locales), @@ -267,7 +267,7 @@ What the `SS_optim` function does: ### Minimum code for running ResistanceGA -The example above outlines how the functions can be used while simulating and generating data with known parameter values. When analyzing your own data, it is not necessary to to use the `r.tran` or `Run_CS` functions. You should only have to use the following functions (although you may want to change settings from defaults). +The example above outlines how the functions can be used while simulating and generating data with known parameter values. When analyzing your own data, it is not necessary to use the `r.tran` or `Run_CS` functions. You should only have to use the following functions (although you may want to change settings from defaults). ```{r eval=FALSE} GA.inputs <- GA.prep(ASCII.dir = write.dir) @@ -326,7 +326,7 @@ t(SS_RESULTS.gdist$ContinuousResults[c(9:11)])) colnames(SS_table) <- c("Truth", "CS.Optimized", "commuteDist.Optimized") SS_table - Truth CS.Optimized commuteDist.Optimized +Truth CS.Optimized commuteDist.Optimized Equation Monomolecular Monomolecular Monomolecular shape 2 1.999999 1.999504 max 275 274.9982 277.0668 @@ -672,7 +672,7 @@ Resist <- Combine_Surfaces(PARM = PARM, Resist$percent.contribution ``` ``` - surface mean l95 u95 +surface mean l95 u95 1 categorical 0.7486327 0.043705035 0.9919852 2 continuous 0.2513673 0.008014751 0.9562950 ``` @@ -748,13 +748,13 @@ Multi.Surface_optim.gd2$percent.contribution ``` ``` > Multi.Surface_optim.gd$percent.contribution - surface mean l95 u95 +surface mean l95 u95 1 categorical 0.7010857 0.020919088 0.9798840 2 continuous 0.2632950 0.017292284 0.9468672 3 feature 0.0356193 0.002216987 0.1878875 > > Multi.Surface_optim.gd2$percent.contribution - surface mean l95 u95 +surface mean l95 u95 1 categorical 0.7174886 0.03268185 0.9841884 2 continuous 0.2825114 0.01581160 0.9673181 ``` @@ -803,14 +803,15 @@ response[lower.tri(response)] <- gd.response ``` ``` - surface avg.R2m avg.weight avg.rank n Percent.top k - -1 categorical.continuous 0.4537753 4.674643e-01 1.436 569 56.9 6 -2 categorical.continuous.feature 0.4545434 4.441112e-01 1.643 404 40.4 8 -3 categorical 0.4374701 8.842401e-02 2.921 27 2.7 4 -4 continuous 0.1533399 1.430606e-07 4.578 0 0.0 4 -5 Distance 0.1512606 1.676945e-07 5.147 0 0.0 2 -6 feature 0.1493708 1.716814e-07 5.275 0 0.0 3 +# A tibble: 6 x 10 + surface avg.AIC avg.AICc avg.weight avg.rank avg.R2m avg.LL n Percent.top k + +1 categorical.continuous 1110 1111 0.468 1.44 0.452 -551 564 56.4 6 +2 cat.cont.feature 1110 1111 0.443 1.63 0.453 -551 410 41.0 8 +3 categorical 1115 1115 0.0892 2.93 0.436 -553 26.0 2.60 4 +4 continuous 1149 1149 0.0000 4.58 0.154 -571 0 0 4 +5 Distance 1150 1150 0.0000 5.16 0.152 -571 0 0 2 +6 feature 1150 1150 0.0000 5.27 0.150 -571 0 0 3 ``` You can see that there is decent support for the composite surface that (incorrectly) includes the feature surface, but that in 57% of bootstrap samples, the composite without the feature surface was the top model. The model selection clearly suggests, correctly, that the multisurface combination with the categorical and continuous surface is the best. There was some ambiguity about whether the categorical.continuous composite was an improvement over the categorical surface alone, but the bootstrap analysis suggests virtually no support for the categorical surface. Further, if you assess the optimized values from this analysis (below), you can see that each surface was parameterized fairly well. The optimized parameters do not perfectly match the true data generating parameters, but overall this resulted in a nearly perfectly correlated surface (r = 1.00). It is worth noting that the optimized surface resulted in a marginal R^2^ = 0.46, suggesting that optimization with 'noisy' genetic data may not result in the highest R^2^ values, but that optimization can nonetheless reliably determine the effects of the landscape on pairwise genetic distance. @@ -822,7 +823,7 @@ row.names(Summary.table)<-c("Category1", "Category2", "Category3", ``` ``` > Summary.table - Truth Optimized +Truth Optimized Category1 1.0 1.00 Category2 250.0 347.56 Category3 100.0 123.80 @@ -849,12 +850,147 @@ pairs(r.stack) * If the optimization ends very quickly (e.g., <40 iterations), you may want to increase the probability of mutation (`pmutation`) and/or the probability of crossover (`pcrossover`). These can be adjusted using `GA.prep`. I have not extensively tested these settings to determine optimal values, but found that the current defaults (`pmutation = 0.125`, `pcrossover = 0.85`) have generally worked quite well with simulated data and produced reproducible estimates with real data. Alternatively, because this is a stochastic optimization, just rerun the optimization (**make sure you have not set a seed!**)\ * Any and all settings of the `ga` function can be adjusted or customized. The main change made from the default setting for optimization of resistance surfaces was to use the "gareal_blxCrossover" method. This greatly improved the search of parameter space.\ * As mentioned above concerning single surface optimization: this is a stochastic optimization process and optimized values will likely differ from run to run. Despite the time involved, it is advised to run all optimizations at least twice to confirm parameter estimates/relative relationship among resistance surfaces. \ -* While there is no established framework for how optimization of resistances surface can or should be done, below is a flowchart of how an analysis might proceed: -![flowchart](figure/FlowChart_Narrow2b.png) \ +# Full Analysis: `all_comb` +There is no established framework for how optimization of resistances surface can or should be done. In dealing with simulated and real data sets, it is often the case that a surface has little to no support when optimized in isolation. But when these 'poorly supported' surfaces are included in multisurface optimizations, they often result in the best fitting model. **As such, I think a step-wise (i.e. single surface followed by selective multi-surface) optimization should be discouraged.** Rather, all possible combinations of surfaces hypothesized to be of importance should be assessed. To facilitate such an analysis, you can use the `all_comb` function. This convenience function will conduct all single surface optimizations, followed by multisurface optimizations up to a specified complexity (i.e. number of surfaces to combine). Following all optimizations, a bootstrap analysis is conducted to assess the level of support for each of the fitted surfaces. The code below will complete the same analysis demonstrated above, but will include all possible combinations of the three resistance surfaces. + +```{r all_comb_demo, eval = F} +# Run `gdist.prep` function as before +gdist.inputs <- gdist.prep(n.Pops = length(sample.locales), + response = gd.response, + samples = sample.locales) + + +# GA.prep for all combinations analysis +# NOTE: Specify "all.comb" for `Results.dir` +GA.inputs <- GA.prep(method = "LL", + ASCII.dir = resistance_surfaces, + Results.dir = "all.comb", + max.cat = 500, + max.cont = 500, + select.trans = list(NA, + "M", # Only monomolecular considered + NA), + parallel = 6) + + +# Replicate optimization runs can completed by specifying +# the a number for `replicate` +all.comb_results <- all_comb(gdist.inputs = gdist.inputs, + GA.inputs = GA.inputs, + results.dir = "C:/ResistanceGA_Examples/all_comb/", + max.combination = 3, + replicate = 1 + ) +``` + +```{r, eval = F} +all.comb_results$summary.table +``` +``` +> all.comb_results$summary.table + Surface obj.func_LL k AIC AICc R2m R2c LL delta.AICc weight +1 categorical -1083.51 4 2175.01 2177.01 0.46 0.46 -1083.51 0.000 0.792 +2 categorical.continuous -1081.73 6 2171.46 2180.12 0.46 0.46 -1081.73 3.109 0.167 +3 categorical.feature -1083.22 6 2174.44 2183.11 0.46 0.46 -1083.22 6.092 0.038 +4 cat.cont.feature -1081.43 8 2170.87 2187.87 0.46 0.46 -1081.43 10.855 0.003 +5 Distance -1114.83 2 2237.66 2234.20 0.16 0.38 -1114.83 57.190 0.000 +6 feature -1114.82 3 2237.63 2236.78 0.15 0.37 -1114.82 59.763 0.000 +7 continuous -1114.47 4 2236.93 2238.93 0.16 0.37 -1114.47 61.919 0.000 +8 continuous.feature -1114.46 5 2236.92 2242.08 0.16 0.37 -1114.46 65.070 0.000 +9 Null -1139.25 1 2284.50 2280.67 0.00 0.20 -1139.25 103.660 0.000 + +``` + + +```{r, eval = F} +all.comb_results$boot.results +``` +``` +> all.comb_results$boot.results +# A tibble: 8 x 10 + surface avg.AIC avg.AICc avg.weight avg.rank avg.R2m n Percent.top k + +1 categorical.continuous 1110 1110 0.33 1.70 0.46 457.0 45.7 6 +2 cat.cont.feature 1109 1110 0.32 2.02 0.46 390.0 39.0 8 +3 categorical 1111 1112 0.18 3.13 0.45 85.0 8.5 4 +4 categorical.feature 1111 1112 0.17 3.15 0.45 68.0 6.8 6 +5 continuous 1149 1149 0.00 5.63 0.15 0.0 0.0 4 +6 continuous.feature 1149 1149 0.00 6.53 0.15 0.0 0.0 5 +7 Distance 1149 1149 0.00 6.88 0.15 0.0 0.0 2 +8 feature 1149 1149 0.00 6.96 0.15 0.0 0.0 3 + +``` +The results from the `all_comb` function are similar to the previous analysis. The frequency that categorical.continuous and categorical.continuous.feature is reduced some because of the addition of categorical.feature, which garnered some support. The preponderance of evidence is still in favor of the categorical.continuous surface as the best resistance surface explaining the observed patterns of genetic differentiation. + +# Fitting MLPE models: `mlpe_rga` +The function `mlpe_rga` can be used to flexibly fit mixed effects models using the MLPE parameterization. In `ResistanceGA`, there is only ever a single predictor variable (pairwise effective distance between sampled populations/individuals). The `mlpe_rga` function allows for the specification of multiple predictor variables. For instance, one may wish to fit pairwise Euclidean distance in addition to effective distance when fitting an MLPE model (e.g., [Row et al. 2017](http://onlinelibrary.wiley.com/doi/10.1002/ece3.2825/abstract?systemMessage=Please+be+advised+that+we+experienced+an+unexpected+issue+that+occurred+on+Saturday+and+Sunday+January+20th+and+21st+that+caused+the+site+to+be+down+for+an+extended+period+of+time+and+affected+the+ability+of+users+to+access+content+on+Wiley+Online+Library.+This+issue+has+now+been+fully+resolved.++We+apologize+for+any+inconvenience+this+may+have+caused+and+are+working+to+ensure+that+we+can+alert+you+immediately+of+any+unplanned+periods+of+downtime+or+disruption+in+the+future. "Row et al.")). +```{r, eval = F} +# Recycling objects from the Example Analysis above: +# Get pairwise Euclidean distances +e.dist <- lower(as.matrix(dist(samples[,2:3]))) + +# Generate to-from object +id <- To.From.ID(nrow(samples)) + +# Create data frame, be sure to add `pop1` from the object +# made by running `To.From.ID` above +df <- data.frame(g.dist = gd.response, + resist.dist = gd.true, + e.dist = e.dist, + pop = id$pop1) + +# Fit/compare model with and without distance +mod1 <- mlpe_rga(formula = g.dist ~ resist.dist + (1 | pop), + data = df) + +mod2 <- mlpe_rga(formula = g.dist ~ resist.dist + e.dist + (1 | pop), + data = df) + +summary(mod2) +AIC(mod1, mod2) + +``` +``` +> summary(mod2) +Linear mixed model fit by maximum likelihood ['lmerMod'] + + AIC BIC logLik deviance df.resid + 2175.4 2193.9 -1082.7 2165.4 295 + +Scaled residuals: + Min 1Q Median 3Q Max +-2.8649 -0.7295 0.0093 0.6209 2.7724 + +Random effects: + Groups Name Variance Std.Dev. + pop (Intercept) 0.00 0.000 + Residual 79.85 8.936 +Number of obs: 300, groups: pop, 25 + +Fixed effects: + Estimate Std. Error t value +(Intercept) 0.48424 2.32249 0.208 +resist.dist 0.99594 0.06830 14.581 +e.dist 0.02051 0.15258 0.134 + +Correlation of Fixed Effects: + (Intr) rsst.d +resist.dist -0.829 +e.dist -0.149 -0.391 + +> AIC(mod1, mod2) + df AIC +mod1 4 2173.429 +mod2 5 2175.411 +``` + +**References** +Row, J. R., S. T. Knick, S. J. Oyler-McCance, S. C. Lougheed, and B. C. Fedy. 2017. Developing approaches for linear mixed modeling in landscape genetics through landscape-directed dispersal simulations. Ecology and Evolution 7:3751--3761. + -### Summary +# Summary Hopefully this vignette/tutorial has demonstrated the functions present in this package and how they can be used together to optimize resistance surfaces in isolation or in combination. These methods require no *a priori* assumptions by the researcher. Optimization is conducted solely on the genetic distance data provided. The goal of this package is to make these methods accessible and useful to others. Development and advancement will continue as long as there is interest and there remains a need. Please contact me () if you encounter issues with any of these functions, need assistance with interpretation, or would like other features added. diff --git a/vignettes/ResistanceGA.html b/vignettes/ResistanceGA.html index 2a5cca7..e1e19f7 100644 --- a/vignettes/ResistanceGA.html +++ b/vignettes/ResistanceGA.html @@ -8,151 +8,72 @@ + + + + +A Vignette/Tutorial to use ResistanceGA - - - - - - - - - - - - - + - - - - -
- - - - - - +

A Vignette/Tutorial to use ResistanceGA

+

Bill Peterman

+

01 February 2018

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
title: “A Vignette/Tutorial to use ResistanceGA”
author: “Bill Peterman”
date: “18 October 2017”
output:
pdf_document: rmarkdown::pdf_document
html_document: rmarkdown::html_vignette
vignette: |
%
%
%

ResistanceGA

@@ -188,7 +109,7 @@

Setup

Install necessary software and packages

If you want to optimize using CIRCUITSCAPE, this package requires that you have CIRCUITSCAPE v4.0 or higher installed on your Windows machine. At this point in time, ResistanceGA can only execute CIRCUITSCAPE on WINDOWS operating systems. You will also need to have R >= v3.0 installed. I would highly recommend installing R studio when working with R.

-

Because it has now been determined that commuteDistance is equivalent to CIRCUITSCAPE, it is highly recommended to to conduct all optimizations using this approach, run in parallel. Doing so allows for all processes to be conducted within R, and will not require execution of software outside of the R environment. If final current flow maps are desired, these can still be generated using the optimized resistance surface using CIRCUITSCAPE.

+

Because it has now been determined that commuteDistance is equivalent to CIRCUITSCAPE, it is highly recommended to conduct all optimizations using this approach, run in parallel. Doing so allows for all processes to be conducted within R, and will not require execution of software outside of the R environment. If final current flow maps are desired, these can still be generated using the optimized resistance surface using CIRCUITSCAPE.

Running on Linux

@@ -203,17 +124,17 @@

Running on Linux

Installation

First, install ResistanceGA from GitHub. This will require the devtools package

-
# Install 'devtools' package, if needed
-if(!("devtools" %in% list.files(.libPaths()))) {
-  install.packages("devtools", repo = "http://cran.rstudio.com", dep = TRUE) 
+
# 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) 

+# Download package, build vignette
+devtools::install_github("wpeterman/ResistanceGA", 
+                         build_vignettes = TRUE)

Load ResistancaGA and clear your workspace.

-
library(ResistanceGA)
-rm(list = ls())
+
library(ResistanceGA)
+rm(list = ls())
@@ -223,16 +144,16 @@

Continuous surface transformations

There are 8 different transformations that can be applied to continuous surfaces. Since the publication of Peterman et al. (2014), I have added Reverse Ricker and Inverse-Reverse Ricker transformation to better cover parameter space. I still think that there are more flexible ways to optimize surfaces, and I’m continuing to develop these as I have time.
Transformations
All of these figures were made with the Plot.trans function. This function returns a ggplot object, which allows you to manipulate some aspects of the plot, as well as determine the resistance value at different levels of your original surface.

-
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))]
+
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))]
## [1] 2.356784
@@ -241,74 +162,74 @@

Example Function Use

Single surface optimization

Make a directory to write ASCII files, CIRCUITSCAPE batch files, and results. It is critical that there are NO SPACES in the specified directory as this will cause functions that interact with CIRCUITSCAPE to fail.

-
if("ResistanceGA_Examples"%in%dir("C:/")==FALSE) 
-  dir.create(file.path("C:/", "ResistanceGA_Examples")) 
+
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")) 
+# 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/"      
+# 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"')

+# 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 resistance surfaces and export as .asc file for use with CIRCUITSCAPE. These surfaces were made using the RandomFields package

-
data(resistance_surfaces)
-continuous <- resistance_surfaces[[2]]
-writeRaster(continuous,
-            filename = paste0(write.dir,"cont.asc"),
-            overwrite = TRUE)
+
data(resistance_surfaces)
+continuous <- resistance_surfaces[[2]]
+writeRaster(continuous,
+            filename = paste0(write.dir,"cont.asc"),
+            overwrite = TRUE)

Load the example sample location data and export as .txt file. This is formatted for input into CIRCUITSCAPE

-
data(samples)
-write.table(samples,file=paste0(write.dir,"samples.txt"),sep="\t",col.names=F,row.names=F)
+
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)])

+# Create a spatial points object for plotting
+sample.locales <- SpatialPoints(samples[,c(2,3)])

Plot surface and overlay the sample points

-
plot(continuous)
-plot(sample.locales, pch = 16, col = "blue", add = TRUE) # Add points
-

+
plot(continuous)
+plot(sample.locales, pch = 16, col = "blue", add = TRUE) # Add points
+

Prepare data for optimization

Run the GA.prep and CS.prep functions

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

Note that response was not defined in CS.prep because it has not been made yet. When doing an actual analysis (not a simulation, as in this example) you will specify your pairwise genetic distance data as the response.

select.trans allows you to specify which transformations can be applied to a surface during optimization. “A” = All; “M” = Monomolecular only; “R” = Ricker only. By default, all transformations will be assessed when a continuous surface is optimized. See GA.prep documentation for more details. In this example, we are constraining the genetic algorithm to only consider Monomolecular transformations. In an actual analysis, such a constraint may or may not be desired.

Transform the raw continuous surface using the Resistance.tran function to apply one of the eight transformations, and then view the transformation using Plot.trans. Note that Plot.trans returns a ggplot2 object as well as the plot. Therefore you can manipulate and modify the plot as desired.

-
r.tran <- Resistance.tran(transformation = "Monomolecular", 
-                          shape = 2, 
-                          max = 275, 
-                          r = continuous) 
-
-plot.t <- Plot.trans(PARM = c(2, 275), 
-                     Resistance = continuous, 
-                     transformation = "Monomolecular")
+
r.tran <- Resistance.tran(transformation = "Monomolecular", 
+                          shape = 2, 
+                          max = 275, 
+                          r = continuous) 
+
+plot.t <- Plot.trans(PARM = c(2, 275), 
+                     Resistance = continuous, 
+                     transformation = "Monomolecular")

Run the transformed resistance surface through CIRCUITSCAPE to get effective resistance between each pair of points. Run.CS returns the lower half of the pairwise resistance matrix for use with the optimization prep functions. This will be our response that we optimize on.

-
# Create the true resistance/response surface
-CS.response <- Run_CS(CS.inputs = CS.inputs,
-                      GA.inputs = GA.inputs,
-                      r = r.tran)
+
# Create the true resistance/response surface
+CS.response <- Run_CS(CS.inputs = CS.inputs,
+                      GA.inputs = GA.inputs,
+                      r = r.tran)

Rerun CS.prep including the newly created CS.response

-
CS.inputs <- CS.prep(n.Pops = length(sample.locales),
-                     response = CS.response,
-                     CS_Point.File = paste0(write.dir,"samples.txt"),
-                     CS.program = CS.program)
+
CS.inputs <- CS.prep(n.Pops = length(sample.locales),
+                     response = CS.response,
+                     CS_Point.File = paste0(write.dir,"samples.txt"),
+                     CS.program = CS.program)

Run the Single surface optimization function (SS_optim). Running this example with the default settings took 147 iterations and ~72 minutes to complete on a computer with an Intel i7 3.6 GHz processor. The data generating values have been precisely recovered.

-
SS_RESULTS <- SS_optim(CS.inputs=CS.inputs,
-                       GA.inputs=GA.inputs)
+
SS_RESULTS <- SS_optim(CS.inputs=CS.inputs,
+                       GA.inputs=GA.inputs)

View the results and compare with truth

SS_table <- data.frame(c("Monomolecular", 2.0, 275),
 t(SS_RESULTS$ContinuousResults[c(3:5)]))
@@ -337,16 +258,16 @@ 
Prepare data for optimization

 * The returned object is a named list containing the tables described above.

 

 
Minimum code for running ResistanceGA

-
The example above outlines how the functions can be used while simulating and generating data with known parameter values. When analyzing your own data, it is not necessary to to use the r.tran or Run_CS functions. You should only have to use the following functions (although you may want to change settings from defaults).

-
GA.inputs <- GA.prep(ASCII.dir = write.dir) 
+
The example above outlines how the functions can be used while simulating and generating data with known parameter values. When analyzing your own data, it is not necessary to use the r.tran or Run_CS functions. You should only have to use the following functions (although you may want to change settings from defaults).

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

+SS_RESULTS <- SS_optim(CS.inputs = CS.inputs,
+                       GA.inputs = GA.inputs)
@@ -357,32 +278,32 @@

Optimization with CIRCUITSCAPE on Linux platform

Optimzation using costDistance (least cost paths) and commuteDistance (equivalent to CS resistance distance)

The above optimization can also be done using functions in gdistance. This approach uses least cost paths or random walk commute times between points. Optimization using gdistance can be substantially faster than optimization with CIRCUITSCAPE. This optimization took 49 iterations and 4 minutes to complete when run in parallel on 4 cores (parallel = 4 in GA.prep). To demonstrate the equivalence of CIRCUITSCAPE and commuteDistance, the simulation below uses the CIRCUITSCAPE resistance distances generated above as the response.

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

Now compare the results from optimization with CIRCUITSCAPE and optimization with commuteDistance

SS_table <- data.frame(c("Monomolecular", 2.0, 275),
 t(SS_RESULTS$ContinuousResults[c(9:11)]),
@@ -396,18 +317,18 @@ 
Optimzation using costDistance (least cost paths) and com
 max                275        274.9982                277.0668

To view the response surface for the Monomolecular optimization of this surface, you can run Grid.Search. This function is only relevant for single continuous surfaces.

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

GRID.Surface  

You can change the color scheme and color breaks by manually recreating the response surface from the generated data [default = topo.colors(20)]

-
filled.contour(Grid.Results$Plot.data,
-               col = rainbow(14),
-               xlab = "Shape parameter",
-               ylab = "Maximum value parameter")
+
filled.contour(Grid.Results$Plot.data,
+               col = rainbow(14),
+               xlab = "Shape parameter",
+               ylab = "Maximum value parameter")

GRID.Surface.update  

Note that actual response surfaces tend to be slightly flatter, and the maximum value for a single surface is more difficult to identify precisely. If you were to add some random noise to the CS.response, the single surface optimization generally would do a good job of recovering the transformation and shape parameters, but the true maximum value may remain elusive. Occasionally the algorithm will get ‘stuck’ trying to optimize on an incorrect transformation. If this happens, rerun the optimization. Of course, you may not know that a surface wasn’t correctly optimized when using real data. For this reason, it is good practice to run all optimizations at least twice to confirm parameter estimates.

@@ -415,85 +336,85 @@

Optimzation using costDistance (least cost paths) and com

Simultaneous optimization of multiple surfaces

First, make a new directory to write ASCII files, CIRCUITSCAPE batch files, and results.

-
if("ResistanceGA_Examples"%in%dir("C:/")==FALSE) 
-  dir.create(file.path("C:/", "ResistanceGA_Examples")) 
+
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")) 
+# 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/"      

+# Directory to write .asc files and results
+write.dir <- "C:/ResistanceGA_Examples/MultipleSurfaces/"

Extract other resistance surfaces from the ‘resistance_surfaces’ raster stack

-
data(resistance_surfaces)
-data(samples)
-sample.locales <- SpatialPoints(samples[ ,c(2,3)])
+
data(resistance_surfaces)
+data(samples)
+sample.locales <- SpatialPoints(samples[ ,c(2,3)])

Visualize each surface:

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

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

Create a raster stack and run the GA.prep function (needed to combine surfaces). In this example, the genetic algorithm will seek to maximize the log likelihood of the fitted model (method = "LL").

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

Transform, reclassify, and combine the three resistance surfaces together. Use an “Inverse-Reverse Monomolecular” transformation of the continuous surface. Visualize this transformation using Plot.trans. The first value of PARM refers to the shape parameter, and the second value refers to the maximum value parameter. Look in the help file for Plot.trans for transformation names/numbers.

-
plot.t <- Plot.trans(PARM = c(3.5, 150),
-                     Resistance = continuous,
-                     transformation = "Inverse-Reverse Monomolecular")
+
plot.t <- Plot.trans(PARM = c(3.5, 150),
+                     Resistance = continuous,
+                     transformation = "Inverse-Reverse Monomolecular")

InverseReverse_Monomolecular
Combine raster surfaces together using Combine_Surfaces. Note that the .asc files are read in alphabetically (if stored in a directory), or else imported in the order they occur in the raster stack. You can check the order of surfaces by inspecting GA.inputs$layer.names. First, define the parameters that will be passed to Combine_Surfaces.

-
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")
+
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")

RESIST.Surface 

Generate new gdist response surface by using Run_gdistance and run gdsit.prep to add response

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

Run MS_optim. Running this multisurface example with the default settings took 215 iterations and ~42 minutes to complete on a computer with an Intel i7 3.6 GHz processor.

-
Multi.Surface_optim <- MS_optim(gdist.inputs = gdist.inputs,
-                                GA.inputs = GA.inputs)
+
Multi.Surface_optim <- MS_optim(gdist.inputs = gdist.inputs,
+                                GA.inputs = GA.inputs)

What the MS_optim function does:
* Read all .asc files that are in the specified ASCII.dir (if using CIRCUITSCAPE), makes a raster stack, and determines whether each is a categorical or continuous surface. A surface is considered categorical if it contains 15 or fewer unique values. * Transformation and resistance values are chosen for each surface, all surfaces are added together, and chosen objective function is obtained from the mixed effects model.
* Several summary outputs are generated
@@ -503,11 +424,11 @@

Simultaneous optimization of multiple surfaces

* In the ‘Plots’ directory there is a 4-panel figure with different model diagnostic plots generated from the fitted mixed effects model of each optimized resistance surface.
* The GA object from the optimization is returned and can be further explored.

The multisurface optimization procedure has done a pretty good job of recovering the relative data generating values. You’ll notice that we have not exactly recovered the values, but that the relative relationship among surfaces is preserved (see below).

-
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")
+
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") 
Summary.table
 Truth Optimized
 Category1        1.0      1.00
@@ -519,44 +440,44 @@ 
Simultaneous optimization of multiple surfaces

 Feature1         1.0      1.00
 Feature2       350.0    326.35

If we rescale both the true and optimized resistance surfaces to have a minimum value of 1, we see that the surfaces are identical. The values for the 3-class categorical surface are the first three values listed, the continuous surface values = 4–6 , and the feature surface values = 7–8. Note that the first value for continuous surfaces identifies the transformation used (the fourth value, here), and is always rounded down (1 = Inverse-Reverse Monomolecular). Visualize and test the equivalence of simulated and optimized resistance surfaces:

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

combined.plots
correlation.plots

If you want to create a CIRCUITSCAPE current surface from either the true or optimized surfaces, this can be done by setting CurrentMap = TRUE and output = "raster" in Run_CS.

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

CS_corr.plot
The optimization can converge on a highly correlated solution, but one that results in relative resistance values that are identical to those of the simulated data. This is important to understand, and interpretation of resistance values should be made with this fact in mind.

@@ -564,64 +485,64 @@

Simultaneous optimization of multiple surfaces

Optimization of a smoothing parameter

Scale is a central theme in landscape ecology, however it seems to be rarely or only indirectly addressed in landscape genetics studies. One way to assess ‘scale’ is to apply a kernel smoothing function to a continuous surface, or to a binary categorical surface. As of version 3.0-4, it is possible to set scale = TRUE in the GA.prep function. Doing so will apply a Gaussian kernel smoothing function to your resistance surface prior to applying a transformation. The sigma parameter is the standard deviation of the Gaussian kernel, and is measured in raster pixels (not map units).

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

smooth.plot
Run optimization to determine the transformation and smoothing parameters.

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

Assess optimized against truth

SS_table <- data.frame(c("Monomolecular", 2.0, 275, 1),
 t(SS_RESULTS.gdist_scale$ContinuousResults[c(9:12)]))
@@ -635,193 +556,193 @@ 
Optimization of a smoothing parameter

 scale                1             0.9916279

Not too bad! As previously mentioned, the maximum value parameter seems to be the hardest to optimize precisely. Note that optimization with a scaling parameter can be done with single surfaces or with a multi-surface composite. However, make sure that only continuous or binary surfaces are included.

-
-

Example Analysis

+
+
+

Example Analysis

Below is an example analysis using the three raster surfaces and sample locations provided with ResistanceGA. The ‘true’ data generating surface in this example will be a combination of the categorical and transformed continuous surface. The feature surface will have no effect on the pairwise distances. This example will show how ResistanceGA can (1) be used for model selection; (2) determine the relative influence of each resistance surface within a composite or multisurface analysis. A full ‘analysis’ is presented below. Note that for the first time in the simulation of data, ‘noise’ is being added to make it more akin to empirical genetic data.

-
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 Reverse Monomolecular transformation of the continuous surface
-# Reverse Monomolecular  = 5
-PARM <- c(1, 250, 100, 1, 1.5, 150)
-# PARM <- c(1, 250, 100, 1, 1.5, 150) # GOOD
-
-
-# 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
-
      surface      mean         l95       u95
+
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

+
surface      mean         l95       u95
 1 categorical 0.7486327 0.043705035 0.9919852
 2  continuous 0.2513673 0.008014751 0.9562950

 
Here, the categorical surface is responsible for ~62% of the total landscape resistance, while the continuous surface is responsible for ~38%.

-
# Create a subdirectory for results
-dir.create(file.path("C:/ResistanceGA_Examples/","run1"))
-dir.create(file.path("C:/ResistanceGA_Examples/","run2"))
+
# 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) 
+# Turn response data into vector
+gd.true <- Run_gdistance(gdist.inputs = gdist.inputs,
+                         r = Resist$combined.surface) 
 
-gd.true <- as.vector(gd.true)
+gd.true <- as.vector(gd.true)
 
-# Add some noise to response
-set.seed(321)
-gd.response <- gd.true + rnorm(length(gd.true), 0, 9)
+# 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
+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
-ecodist::mantel(gd.response ~ gd.true + as.vector(dist(samples[, c(2, 3)]))) # Partial  Mantel r = 0.64
-plot(gd.response ~ as.vector(dist(samples[, c(2, 3)])), xlab = "Euclidean distance")

-
Noisy Genetic x Resistance distance

-Resistance dist x Euc. dist

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

Noisy Genetic x Resistance distance
+Resistance dist x Euc. dist

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

What are the Percent contributions of individual surfaces to the combined surface?

-
Multi.Surface_optim.gd$percent.contribution
+
Multi.Surface_optim.gd$percent.contribution
 
-Multi.Surface_optim.gd2$percent.contribution

+Multi.Surface_optim.gd2$percent.contribution
> Multi.Surface_optim.gd$percent.contribution
-      surface      mean         l95       u95
+surface      mean         l95       u95
 1 categorical 0.7010857 0.020919088 0.9798840
 2  continuous 0.2632950 0.017292284 0.9468672
 3     feature 0.0356193 0.002216987 0.1878875
 > 
 > Multi.Surface_optim.gd2$percent.contribution
-      surface      mean        l95       u95
+surface      mean        l95       u95
 1 categorical 0.7174886 0.03268185 0.9841884
 2  continuous 0.2825114 0.01581160 0.9673181

We can see that the feature surface has a minimal contribution to the total resistance, and that the mean contributions of the categorical and continuous surfaces mirror the ‘true’ contributions (categorical = 75%, continuous = 25%).

-

Bootstrap Analysis

It is often the case that there is weak or ambiguous support for a single best resistance surface (i.e. similar AIC scores), as can be seen in the AIC table below. Based strictly on AIC or marginal R2, there really is no difference between the multisurface optimization with and without the feature surface included. If we account for the additional parameters (k) present when including the feature surface, AICc does pretty clearly suggest that the composite surface without the feature is best supported.

-
Surface                         obj.func_LL   k   AIC         AICc      R2m   R2c     LL
-categorical.continuous          -1081.73        6     2171.46     2180.12     0.46  0.46    -1081.73
-categorical                     -1086.49        4     2180.98     2182.98     0.44  0.44    -1086.49
-categorical.continuous.feature  -1081.44        8     2170.88     2187.88     0.46  0.46    -1081.44
-Distance                          -1114.83      2     2237.66     2234.20     0.16  0.38    -1114.83
-feature                         -1114.82        3     2237.63     2236.78     0.15  0.37    -1114.82
-continuous                      -1114.47        4     2236.93     2238.93     0.16  0.37    -1114.47
-Null                              -1139.25      1     2284.50     2280.67     0.00  0.20    -1139.25
+
Surface                       obj.func_LL   k     AIC         AICc      R2m   R2c     LL
+categorical.continuous          -1081.73  6   2171.46     2180.12     0.46  0.46    -1081.73
+categorical                     -1086.49    4     2180.98     2182.98     0.44  0.44    -1086.49
+categorical.continuous.feature  -1081.44    8     2170.88     2187.88     0.46  0.46    -1081.44
+Distance                          -1114.83  2     2237.66     2234.20     0.16  0.38    -1114.83
+feature                         -1114.82    3     2237.63     2236.78     0.15  0.37    -1114.82
+continuous                      -1114.47    4     2236.93     2238.93     0.16  0.37    -1114.47
+Null                              -1139.25  1     2284.50     2280.67     0.00  0.20    -1139.25

Another approach that may lead to additional insight beyond the AIC(c) values generated during optimization is a bootstrap analysis using the Resist.boot function. This approach will assess the relative support for each optimized resistance surface through a ’pseudo’bootstrap where the sample locations/individuals as well as the resistance distance matrices are sub-sampled at each bootstrap iteration (without replacement), and the MLPE model is refit and the AIC scores calculated. Resistance surfaces are NOT re-optimized during this process. This analysis will assess the robustness of the optimized resistance surface given different combinations of samples. If the observed patterns in the resistance-genetics relationship are driven by one or a few sample locations, this analysis may reveal this.

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

-
                         surface   avg.R2m   avg.weight avg.rank     n Percent.top     k
-                           <chr>     <dbl>        <dbl>    <dbl> <dbl>       <dbl> <dbl>
-1         categorical.continuous 0.4537753 4.674643e-01    1.436   569        56.9     6
-2 categorical.continuous.feature 0.4545434 4.441112e-01    1.643   404        40.4     8
-3                    categorical 0.4374701 8.842401e-02    2.921    27         2.7     4
-4                     continuous 0.1533399 1.430606e-07    4.578     0         0.0     4
-5                       Distance 0.1512606 1.676945e-07    5.147     0         0.0     2
-6                        feature 0.1493708 1.716814e-07    5.275     0         0.0     3

+)
+
# A tibble: 6 x 10
+  surface                  avg.AIC avg.AICc  avg.weight avg.rank avg.R2m avg.LL   n     Percent.top  k
+  <fctr>                    <dbl>    <dbl>    <dbl>     <dbl>    <dbl>   <dbl>   <dbl>   <dbl>      <dbl>
+1 categorical.continuous    1110     1111    0.468      1.44    0.452    -551    564     56.4       6
+2 cat.cont.feature          1110     1111    0.443      1.63    0.453    -551    410     41.0       8
+3 categorical               1115     1115    0.0892     2.93    0.436    -553    26.0     2.60      4
+4 continuous                1149     1149    0.0000     4.58    0.154    -571     0       0         4
+5 Distance                  1150     1150    0.0000     5.16    0.152    -571     0       0         2
+6 feature                   1150     1150    0.0000     5.27    0.150    -571     0       0         3

You can see that there is decent support for the composite surface that (incorrectly) includes the feature surface, but that in 57% of bootstrap samples, the composite without the feature surface was the top model.
The model selection clearly suggests, correctly, that the multisurface combination with the categorical and continuous surface is the best. There was some ambiguity about whether the categorical.continuous composite was an improvement over the categorical surface alone, but the bootstrap analysis suggests virtually no support for the categorical surface. Further, if you assess the optimized values from this analysis (below), you can see that each surface was parameterized fairly well. The optimized parameters do not perfectly match the true data generating parameters, but overall this resulted in a nearly perfectly correlated surface (r = 1.00). It is worth noting that the optimized surface resulted in a marginal R2 = 0.46, suggesting that optimization with ‘noisy’ genetic data may not result in the highest R2 values, but that optimization can nonetheless reliably determine the effects of the landscape on pairwise genetic distance.

-
Summary.table <- data.frame(PARM,round(t(Multi.Surface_optim.gd2$GA.summary@solution),2))
-colnames(Summary.table)<-c("Truth", "Optimized")
-row.names(Summary.table)<-c("Category1", "Category2", "Category3", 
-                            "Transformation", "Shape", "Max")
+
Summary.table <- data.frame(PARM,round(t(Multi.Surface_optim.gd2$GA.summary@solution),2))
+colnames(Summary.table)<-c("Truth", "Optimized")
+row.names(Summary.table)<-c("Category1", "Category2", "Category3", 
+                            "Transformation", "Shape", "Max") 
> Summary.table
-               Truth Optimized
+Truth Optimized
 Category1        1.0      1.00
 Category2      250.0    347.56
 Category3      100.0    123.80
 Transformation   1.0      1.60
 Shape            1.5      3.03
 Max            150.0     83.75
-
opt.r <- raster("C:/ResistanceGA_Examples/run2/Results/categorical.continuous.asc")
-r.stack <- stack(Resist$combined.surface, opt.r) 
-names(r.stack) <- c("Truth", "Optimized")
+
opt.r <- raster("C:/ResistanceGA_Examples/run2/Results/categorical.continuous.asc")
+r.stack <- stack(Resist$combined.surface, opt.r) 
+names(r.stack) <- c("Truth", "Optimized")
 
-plot(r.stack)
+plot(r.stack)
 
-pairs(r.stack)

+pairs(r.stack)

Side by Side
Example Correlation

@@ -830,37 +751,134 @@

Bootstrap Analysis

* If the optimization ends very quickly (e.g., <40 iterations), you may want to increase the probability of mutation (pmutation) and/or the probability of crossover (pcrossover). These can be adjusted using GA.prep. I have not extensively tested these settings to determine optimal values, but found that the current defaults (pmutation = 0.125, pcrossover = 0.85) have generally worked quite well with simulated data and produced reproducible estimates with real data. Alternatively, because this is a stochastic optimization, just rerun the optimization (make sure you have not set a seed!)
* Any and all settings of the ga function can be adjusted or customized. The main change made from the default setting for optimization of resistance surfaces was to use the “gareal_blxCrossover” method. This greatly improved the search of parameter space.
* As mentioned above concerning single surface optimization: this is a stochastic optimization process and optimized values will likely differ from run to run. Despite the time involved, it is advised to run all optimizations at least twice to confirm parameter estimates/relative relationship among resistance surfaces.
-* While there is no established framework for how optimization of resistances surface can or should be done, below is a flowchart of how an analysis might proceed:

-

flowchart
-

-
-

Summary

-

Hopefully this vignette/tutorial has demonstrated the functions present in this package and how they can be used together to optimize resistance surfaces in isolation or in combination. These methods require no a priori assumptions by the researcher. Optimization is conducted solely on the genetic distance data provided. The goal of this package is to make these methods accessible and useful to others. Development and advancement will continue as long as there is interest and there remains a need. Please contact me (bill.peterman@gmail.com) if you encounter issues with any of these functions, need assistance with interpretation, or would like other features added.

+# Full Analysis: all_comb There is no established framework for how optimization of resistances surface can or should be done. In dealing with simulated and real data sets, it is often the case that a surface has little to no support when optimized in isolation. But when these ‘poorly supported’ surfaces are included in multisurface optimizations, they often result in the best fitting model. As such, I think a step-wise (i.e. single surface followed by selective multi-surface) optimization should be discouraged. Rather, all possible combinations of surfaces hypothesized to be of importance should be assessed. To facilitate such an analysis, you can use the all_comb function. This convenience function will conduct all single surface optimizations, followed by multisurface optimizations up to a specified complexity (i.e. number of surfaces to combine). Following all optimizations, a bootstrap analysis is conducted to assess the level of support for each of the fitted surfaces. The code below will complete the same analysis demonstrated above, but will include all possible combinations of the three resistance surfaces.

+
# Run `gdist.prep` function as before
+gdist.inputs <- gdist.prep(n.Pops = length(sample.locales),
+                           response = gd.response,
+                           samples = sample.locales)
+
+
+# GA.prep for all combinations analysis
+# NOTE: Specify "all.comb" for `Results.dir`
+GA.inputs <- GA.prep(method = "LL",
+                     ASCII.dir = resistance_surfaces,
+                     Results.dir =  "all.comb", 
+                     max.cat = 500,
+                     max.cont = 500,
+                     select.trans = list(NA, 
+                                         "M", # Only monomolecular considered
+                                         NA),
+                     parallel = 6)
+
+
+# Replicate optimization runs can completed by specifying
+# the a number for `replicate`
+all.comb_results <- all_comb(gdist.inputs = gdist.inputs,
+                             GA.inputs = GA.inputs,
+                             results.dir = "C:/ResistanceGA_Examples/all_comb/",
+                             max.combination = 3,
+                             replicate = 1
+                             )
+
all.comb_results$summary.table
+
> all.comb_results$summary.table
+    Surface              obj.func_LL k  AIC     AICc     R2m  R2c       LL    delta.AICc weight
+1   categorical             -1083.51  4 2175.01 2177.01 0.46  0.46  -1083.51    0.000    0.792
+2   categorical.continuous  -1081.73  6 2171.46 2180.12 0.46  0.46  -1081.73    3.109    0.167
+3   categorical.feature     -1083.22  6 2174.44 2183.11 0.46  0.46  -1083.22    6.092    0.038
+4   cat.cont.feature        -1081.43  8 2170.87 2187.87 0.46  0.46  -1081.43    10.855   0.003
+5   Distance                -1114.83  2 2237.66 2234.20 0.16  0.38  -1114.83    57.190   0.000
+6   feature                 -1114.82  3 2237.63 2236.78 0.15  0.37  -1114.82    59.763   0.000
+7   continuous              -1114.47  4 2236.93 2238.93 0.16  0.37  -1114.47    61.919   0.000
+8   continuous.feature      -1114.46  5 2236.92 2242.08 0.16  0.37  -1114.46    65.070   0.000
+9   Null                    -1139.25  1 2284.50 2280.67 0.00  0.20  -1139.25    103.660  0.000
+
+
all.comb_results$boot.results
+
> all.comb_results$boot.results
+# A tibble: 8 x 10
+                 surface    avg.AIC avg.AICc avg.weight avg.rank    avg.R2m   n     Percent.top k
+                  <chr>    <dbl>     <dbl>   <dbl>         <dbl>      <dbl>  <dbl>    <dbl>     <dbl>
+1   categorical.continuous  1110    1110     0.33           1.70      0.46   457.0    45.7      6
+2   cat.cont.feature        1109    1110     0.32           2.02      0.46   390.0    39.0      8
+3   categorical             1111    1112     0.18           3.13      0.45   85.0      8.5      4
+4   categorical.feature     1111    1112     0.17           3.15      0.45   68.0      6.8      6
+5   continuous              1149    1149     0.00           5.63      0.15   0.0       0.0      4
+6   continuous.feature      1149    1149     0.00           6.53      0.15   0.0       0.0      5
+7   Distance                1149    1149     0.00           6.88      0.15   0.0       0.0      2
+8   feature                 1149    1149     0.00           6.96      0.15   0.0       0.0      3
+
+

The results from the all_comb function are similar to the previous analysis. The frequency that categorical.continuous and categorical.continuous.feature is reduced some because of the addition of categorical.feature, which garnered some support. The preponderance of evidence is still in favor of the categorical.continuous surface as the best resistance surface explaining the observed patterns of genetic differentiation.

+
+
+

Fitting MLPE models: mlpe_rga

+

The function mlpe_rga can be used to flexibly fit mixed effects models using the MLPE parameterization. In ResistanceGA, there is only ever a single predictor variable (pairwise effective distance between sampled populations/individuals). The mlpe_rga function allows for the specification of multiple predictor variables. For instance, one may wish to fit pairwise Euclidean distance in addition to effective distance when fitting an MLPE model (e.g., Row et al. 2017).

+
# Recycling objects from the Example Analysis above:
+# Get pairwise Euclidean distances
+e.dist <- lower(as.matrix(dist(samples[,2:3])))
+
+# Generate to-from object
+id <- To.From.ID(nrow(samples))
+
+# Create data frame, be sure to add `pop1` from the object
+# made by running `To.From.ID` above
+df <- data.frame(g.dist = gd.response,
+                 resist.dist = gd.true,
+                 e.dist = e.dist,
+                 pop = id$pop1)
+
+# Fit/compare model with and without distance
+mod1 <- mlpe_rga(formula = g.dist ~ resist.dist + (1 | pop),
+                 data = df)
+
+mod2 <- mlpe_rga(formula = g.dist ~ resist.dist + e.dist + (1 | pop),
+                 data = df)
+
+summary(mod2)
+AIC(mod1, mod2)
+
> summary(mod2)
+Linear mixed model fit by maximum likelihood  ['lmerMod']
+
+     AIC      BIC   logLik deviance df.resid 
+  2175.4   2193.9  -1082.7   2165.4      295 
+
+Scaled residuals: 
+    Min      1Q  Median      3Q     Max 
+-2.8649 -0.7295  0.0093  0.6209  2.7724 
+
+Random effects:
+ Groups   Name        Variance Std.Dev.
+ pop      (Intercept)  0.00    0.000   
+ Residual             79.85    8.936   
+Number of obs: 300, groups:  pop, 25
+
+Fixed effects:
+            Estimate Std. Error t value
+(Intercept)  0.48424    2.32249   0.208
+resist.dist  0.99594    0.06830  14.581
+e.dist       0.02051    0.15258   0.134
+
+Correlation of Fixed Effects:
+            (Intr) rsst.d
+resist.dist -0.829       
+e.dist      -0.149 -0.391
+
+> AIC(mod1, mod2)
+     df      AIC
+mod1  4 2173.429
+mod2  5 2175.411
+

References
+Row, J. R., S. T. Knick, S. J. Oyler-McCance, S. C. Lougheed, and B. C. Fedy. 2017. Developing approaches for linear mixed modeling in landscape genetics through landscape-directed dispersal simulations. Ecology and Evolution 7:3751–3761.

+
+
+

Summary

+

Hopefully this vignette/tutorial has demonstrated the functions present in this package and how they can be used together to optimize resistance surfaces in isolation or in combination. These methods require no a priori assumptions by the researcher. Optimization is conducted solely on the genetic distance data provided. The goal of this package is to make these methods accessible and useful to others. Development and advancement will continue as long as there is interest and there remains a need. Please contact me (bill.peterman@gmail.com) if you encounter issues with any of these functions, need assistance with interpretation, or would like other features added.

Acknowledgements

A huge thanks to Grant Connette for many discussions related to development and implementation of these methods!

-
- - - - - -