Skip to content

Commit

Permalink
changed com_data and added Distance matrix option for fit_models
Browse files Browse the repository at this point in the history
  • Loading branch information
derek-corcoran-barrios committed Jul 4, 2023
1 parent 16c0979 commit 9a379d1
Show file tree
Hide file tree
Showing 4 changed files with 66 additions and 32 deletions.
4 changes: 3 additions & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ Description: Provides tools for model selection and model averaging of PerMANOVA
(<doi:10.1007/b97636>).
License: MIT + file LICENSE
Encoding: UTF-8
Roxygen: list(markdown = TRUE)
Roxygen: list(markdown = TRUE, roclets = c ("namespace", "rd", "srr::srr_stats_roclet"))
RoxygenNote: 7.2.3
Imports:
broom,
Expand All @@ -41,3 +41,5 @@ Suggests:
testthat (>= 3.0.0)
Config/testthat/edition: 3
VignetteBuilder: knitr
Enhances:
betapart
68 changes: 43 additions & 25 deletions R/fit_models.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,17 @@
#' @description This function fits PERMANOVA models for all combinations of variables in a given dataset, and arranges the models by Akaike Information Criterion (AICc) score. The function also calculates the maximum variance inflation factor (max_vif) for each model.
#'
#' @param all_forms A data frame generated by \code{\link{make_models}}
#' @param veg_data A dataset with vegetation presence absense or abundance data
#' @param com_data A dataset with community presence absense or
#' abundance data, you can also use a dist class file generated from
#' vegan, betapart or other packages
#' @param env_data A dataset with the variables described in all_froms
#' @param ncores An integer specifying the number of cores to use for parallel processing
#' @param log logical if true, a log file will be generated
#' @param verbose logical, defaults TRUE, sends messages about processing times
#' @param logfile the text file that will be generated as a log
#' @param multiple after how many loops to write a log file
#' @param method method for distance from \code{\link{vegdist}}
#' @param method method for distance from \code{\link{vegdist}},
#' this will be ignored if com_data is a distance file
#' @param strata a block variable similar to the use in \code{\link{adonis2}}
#'
#' @return A data.frame with fitted models arranged by AICc, including the formula used, the number of
Expand All @@ -28,6 +31,7 @@
#' @importFrom stats rnorm lm as.formula complete.cases
#' @examples
#' \donttest{
#' ## example with dataframe as community data
#' library(vegan)
#' data(dune)
#' data(dune.env)
Expand All @@ -36,9 +40,19 @@
#'
#' fit_models(
#' all_forms = AllModels,
#' veg_data = dune,
#' com_data = dune,
#' env_data = dune.env
#' )
#' ## example with distance as community data
#' library(betapart)
#' Distance <- beta.pair.abund(dune)
#' Distance <- Distance$beta.bray.bal
#' fit_models(
#' all_forms = AllModels,
#' com_data = Distance,
#' env_data = dune.env
#' )
#'
#' }
#'
#' @references
Expand All @@ -49,7 +63,7 @@
#' @export

fit_models <- function(all_forms,
veg_data,
com_data,
env_data,
method = "bray",
ncores = 2,
Expand All @@ -69,32 +83,38 @@ fit_models <- function(all_forms,
if (!("max_vif" %in% colnames(meta_data))) {
meta_data$max_vif <- NA
}
vegetation_data <- veg_data

# Check for missing values
missing_rows <- !complete.cases(env_data)

if (any(missing_rows)) {
if (verbose) {
# Print message about missing rows and columns
message(sprintf("Removing %d rows with missing values\n", sum(missing_rows)))
message("Columns with missing values: ")
message(names(env_data)[colSums(is.na(env_data)) > 0], sep = ", ")
community_data <- com_data

if(class(community_data) != "dist"){
# Check for missing values
missing_rows <- !complete.cases(env_data)

if (any(missing_rows)) {
if (verbose) {
# Print message about missing rows and columns
message(sprintf("Removing %d rows with missing values\n", sum(missing_rows)))
message("Columns with missing values: ")
message(names(env_data)[colSums(is.na(env_data)) > 0], sep = ", ")
}
}
}

# Filter out missing rows
new_env_data <- env_data[complete.cases(env_data), ]
vegetation_data <- vegetation_data[complete.cases(env_data), ]
# Filter out missing rows
env_data <- env_data[complete.cases(env_data), ]
community_data <- community_data[complete.cases(env_data), ]
Distance <- vegan::vegdist(community_data, method = method)
}

if(class(community_data) != "dist"){
Distance <- community_data
}

cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)

Distance <- vegan::vegdist(vegetation_data, method = method)

Fs <- foreach(x = 1:nrow(meta_data), .packages = c("vegan", "dplyr", "AICcPermanova", "tidyr", "broom"), .combine = bind_rows, .export = c("Distance")) %dopar% {
Response <- new_env_data

Response <- env_data
Response$y <- rnorm(n = nrow(Response))
gc()

Expand All @@ -110,8 +130,6 @@ fit_models <- function(all_forms,
Model <- try(with(Response, vegan::adonis2(as.formula(Temp$form[1]), data = Response, by = "margin", strata = strata_factor)), silent = TRUE)
}



Temp <- tryCatch(
expr = cbind(Temp, AICcPermanova::AICc_permanova2(Model)),
error = function(e) NA
Expand Down Expand Up @@ -142,7 +160,6 @@ fit_models <- function(all_forms,
}
)


if (log) {
if ((x %% multiple) == 0) {
sink(logfile, append = TRUE)
Expand All @@ -152,7 +169,8 @@ fit_models <- function(all_forms,
}
}

Temp <- bind_cols(Temp, Rs)
Temp <- dplyr::bind_cols(Temp, Rs)
sink(logfile, append = TRUE)
Temp
}
parallel::stopCluster(cl)
Expand Down
22 changes: 18 additions & 4 deletions man/fit_models.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

4 changes: 2 additions & 2 deletions vignettes/getting-started-with-aiccpermanova.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ After filtering out collinear models, we can fit all the remaining non-collinear
```{r}
Fitted <- fit_models(
all_forms = NonColinear,
veg_data = dune,
com_data = dune,
env_data = dune.env,
ncores = 1,
method = "bray"
Expand All @@ -134,7 +134,7 @@ If there is a block variable to be used (such as Use in the dune.env object), yo
```{r fitwihtstrata}
Fitted2 <- fit_models(
all_forms = NonColinear,
veg_data = dune,
com_data = dune,
env_data = dune.env,
ncores = 1,
method = "bray",
Expand Down

0 comments on commit 9a379d1

Please sign in to comment.