Merge pull request #5 from wpeterman/dev_3.1-4

Merge of development version 3.1-4

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: ResistanceGA
 Title: Optimize resistance surfaces using genetic algorithms
-Version: 3.1-3
+Version: 4.0
 Author: Bill Peterman <[email protected]>
 Maintainer: Bill Peterman <[email protected]>
 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

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,16 +54,22 @@ 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)
 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)

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.

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 <Bill.Peterman@@gmail.com>
-#' @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)
 }

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 <Bill.Peterman@@gmail.com>
 Combine_Surfaces <-

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 <Bill.Peterman@@gmail.com>
 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],