R/AICc_ens_metrics.R

In collinmmccabe/enss: Effective Network Size Simulation

#' AICc_ens_metrics
#'
#' AICc_ens_metrics calculates the best fit model out of all possible linear combinations of metrics as explanatory variables and ENS as the response.
#'
#' @param ens a vector of integer values indicating the effective network size of each observed network, generated by iteratively running the enss function estimate_ens() or estimate_backtrans_ens() for each edgelist generated by the function import_emp().
#' @param metrics a dataframe generated by the enss function calculate_metrics() containing metri The dataframe of metrics should pair with the vector of ens values, such that metrics[i, ] is equal to ens[i] and the number of rows in the metrics dataframe is equal to the length of the ens vector.
#'
#' @return either a single model of best fit, with coefficients for each variable included in the model or a list of models and their coefficients, if there was a tie among AICc values.
#' @export
#'
#' @examples
#' emp_el <- readRDS("./data/Effective_el_emp")
#' emp_el <- list(emp_el[[1]], emp_el[[3]], emp_el[[8]], emp_el[[9]], emp_el[[10]], emp_el[[12]], emp_el[[14]], emp_el[[16]], emp_el[[20]], emp_el[[21]], emp_el[[22]], emp_el[[23]], emp_el[[30]], emp_el[[26]], emp_el[[29]], emp_el[[13]], emp_el[[34]], emp_el[[38]], emp_el[[39]], emp_el[[41]])
#' unweighted_metrics <- calculate_metrics(emp_el)
#' ens_SI <- c(7, 36, 20, 9, 11, 43, 13, 16, 16, 31, 36, 20, 34, 25, 11, 10, 24, 10, 9, 15)
#' SI_AICc <- AICc_ens_metrics(ens_SI, unweighted_metrics)
AICc_ens_metrics <- function(ens, metrics) {

  # Code adapted from ryouready.wordpress.com/2009/02/06/r-calculating-all-possible-linear-regression-models-for-a-given-set-of-predictors/

  regressors <- colnames(metrics[-1])

  data <- cbind(metrics, ens = ens)

  regMat <- expand.grid(c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE))
  regMat <- regMat[-(dim(regMat)[1]),]

  names(regMat) <- regressors

  allModelsList <- apply(regMat, 1, function(x) as.formula(paste(c("ens ~ ", paste(regressors[x], collapse=" + ")), collapse="")))

  allModelsResults <- lapply(allModelsList, function(x) lm(x, data=data))

  AICc_vector <- c(NA)
  for(i in 1:length(allModelsResults)) {

    AICc_vector[i] <- MuMIn::AICc(allModelsResults[[i]])
  }

  AICc_returns <- list(NA); k <- 1
  for(j in which(AICc_vector < AICc_vector[which.min(AICc_vector)] + 2)) {

    AICc_returns[[k]] <- allModelsResults[[j]]
    k <- k + 1
  }

  return(AICc_returns)
}