From 98b5827c9049da53315c7031873527636e1e31ba Mon Sep 17 00:00:00 2001 From: "Adam B. Smith" Date: Thu, 5 Dec 2024 16:32:00 -0600 Subject: [PATCH 1/6] Allow users to remove unconverged/boundary models --- R/trainGlm.r | 120 +++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 108 insertions(+), 12 deletions(-) diff --git a/R/trainGlm.r b/R/trainGlm.r index 2875601..1ef27fd 100644 --- a/R/trainGlm.r +++ b/R/trainGlm.r @@ -25,6 +25,8 @@ #' \item A numeric vector of weights, one per row in \code{data}. #' \item The name of the column in \code{data} that contains site weights. #' } +#' @param removeInvalid Logical. If \code{TRUE} (default), remove models that either did not converge or have parameter estimates near the boundaries (usually negative or positive infinity). If you run this function with `construct = TRUE` (i.e., construct a "full" model from the best "small" models), then any small model that either did not converge or had parameters that are near the boundary (usually negative or positive infinity) are removed from consideration as terms in "full" model. +#' @param failIfInvalid Logical. If \code{TRUE} (default), and the "full" model either does not converge or has parameters near the boundary, then the function will fail. If \code{FALSE}, then return \code{NULL} in this case. #' @param out Character vector. One or more values: #' \itemize{ #' \item \code{'model'}: Model with the lowest AICc. @@ -35,7 +37,7 @@ #' @param verbose Logical. If \code{TRUE} then display progress. #' @param ... Arguments to pass to \code{glm}. #' -#' @returns The object that is returned depends on the value of the \code{out} argument. It can be a model object, a data frame, a list of models, or a list of all two or more of these. If \code{scale} is \code{TRUE}, any model object will also have an element named \code{$scale}, which contains the means and standard deviations for predictors that are not factors. +#' @returns The object that is returned depends on the value of the \code{out} argument. It can be a model object, a data frame, a list of models, or a list of all two or more of these. If \code{scale} is \code{TRUE}, any model object will also have an element named \code{$scale}, which contains the means and standard deviations for predictors that are not factors. The data frame reports the AICc for all of the models, sorted by best to worst. The \code{converged} column indicates whether the model converged ("\code{TRUE}" is good), and the \code{boundary} column whether the model parameters are near the boundary (usually, negative or positive infinity; "\code{FALSE}" is good). #' #' @seealso \code{\link[stats]{glm}} #' @@ -58,6 +60,8 @@ trainGLM <- function( maxTerms = 8, w = TRUE, family = stats::binomial(), + removeInvalid = TRUE, + failIfInvalid = TRUE, out = 'model', cores = 1, verbose = FALSE, @@ -79,6 +83,8 @@ trainGLM <- function( presPerTermFinal <- 10 maxTerms <- 8 w <- TRUE + removeInvalid <- TRUE + failIfInvalid <- TRUE scale <- TRUE family <- stats::binomial() out <- 'model' @@ -191,11 +197,29 @@ trainGLM <- function( ) } - assess <- assess[order(assess$AICc), , drop=FALSE] + if (removeInvalid) { + + bads <- which(!assess$converged | assess$boundary) + if (length(bads) > 0) { + assess <- assess[-bads, , drop = FALSE] + + if (nrow(assess) == 0) { + msg <- 'No single-term models converged or all models had parameter estimates near the boundary.' + if (failIfInvalid) { + stop(msg) + } else { + warning(msg) + return(NULL) + } + } + } + } + + assess <- assess[order(assess$AICc), , drop = FALSE] rownames(assess) <- NULL if (verbose) { - omnibus::say('Term-by-term evaluation:', pre=1) + omnibus::say('Term-by-term evaluation:', pre = 1) print(assess) utils::flush.console() } @@ -232,7 +256,7 @@ trainGLM <- function( numTerms <- lengths(regmatches(thisForm, gregexpr('\\+', thisForm))) + 1 start <- rep(0, numTerms) - model <- stats::glm( + model <- suppressWarnings(stats::glm( formula = stats::as.formula(thisForm), family = family, data = data, @@ -240,23 +264,43 @@ trainGLM <- function( weights = w, start = start, ... - ) + )) AICc <- AICcmodavg::AICc(model) tuning <- data.frame( model = form, + converged = model$converged, + boundary = model$boundary, AICc = AICc ) models <- NULL + if (!tuning$converged | tuning$boundary) { + + msg <- 'The model did not converge and/or estimates are near boundary conditions.' + if (removeInvalid) { + + if (failIfInvalid) { + stop(msg) + } else { + warning(msg) + return(NULL) + } + } else { + warning(msg) + } + + } + + if (verbose) { omnibus::say('Final model (construction from best terms, but no selection):', pre=1) print(summary(model)) utils::flush.console() - + } ### model selection @@ -373,15 +417,44 @@ trainGLM <- function( } # compile all possible models and rank - tuning <- data.frame(model = rep(NA_character_, length(assess)), AICc = rep(NA_real_, length(assess))) + n <- length(assess) + tuning <- data.frame( + model = rep(NA_character_, n), + converged = rep(NA, n), + boundary = rep(NA, n), + AICc = rep(NA_real_, n) + ) + models <- list() for (i in seq_along(assess)) { models[[i]] <- assess[[i]]$model if (!is.na(scale)) if (scale) models[[i]]$scale <- scales tuning$model[i] <- assess[[i]]$formula + tuning$converged[i] <- assess[[i]]$converged + tuning$boundary[i] <- assess[[i]]$boundary tuning$AICc[i] <- assess[[i]]$AICc } + if (removeInvalid) { + + bads <- which(!tuning$converged | tuning$boundary) + if (length(bads) > 0) { + + tuning <- tuning[-bads, , drop = FALSE] + models <- models[-bads] + + if (nrow(tuning) == 0) { + msg <- 'No models converged or all had parameter estimates near the boundary of parameter space.' + if (failIfInvalid) { + stop(msg) + } else { + warning(msg) + return(NULL) + } + } + } + } + ranks <- order(tuning$AICc) models <- models[ranks] tuning <- tuning[ranks, , drop=FALSE] @@ -409,7 +482,7 @@ trainGLM <- function( numTerms <- lengths(regmatches(thisForm, gregexpr('\\+', thisForm))) + 1 start <- rep(0, numTerms) - model <- stats::glm( + model <- suppressWarnings(stats::glm( formula = stats::as.formula(thisForm), family = family, data = data, @@ -417,19 +490,37 @@ trainGLM <- function( weights = w, start = start, ... - ) + )) AICc <- AICcmodavg::AICc(model) tuning <- data.frame( model = form, + converged = model$converged, + boundary = model$boundary, AICc = AICc ) models <- NULL - if (select) warning('Model selection is not performed when argument "construct" is FALSE.') + if (select) warning('Model selection is not performed when argument `construct` is FALSE.') + if (!model$converged | model$boundary) { + + msg <- 'The model did not converge and/or parameters are near the boundary space.' + if (removeInvalid) { + + if (failIfInvalid) { + stop(msg) + } else { + warning(msg) + return(NULL) + } + } else { + warning(msg) + } + } + if (verbose) { omnibus::say('Model (no construction or selection):', level=2) @@ -454,6 +545,7 @@ trainGLM <- function( output <- models } else if ('model' %in% out) { output <- model + if (!tuning$converged[1] | tuning$boundary[1]) warning('The top model did not converge and/or had parameter estimates near the boundary.') } else if ('tuning' %in% out) { output <- tuning } @@ -500,7 +592,7 @@ trainGLM <- function( numTerms <- lengths(regmatches(thisForm, gregexpr('\\+', thisForm))) + 1 start <- rep(0, numTerms) - model <- stats::glm( + model <- suppressWarnings(stats::glm( formula = stats::as.formula(thisForm), data = data, family = family, @@ -508,7 +600,7 @@ trainGLM <- function( method = method, start = start, ... - ) + )) AICc <- AICcmodavg::AICc(model) @@ -519,6 +611,8 @@ trainGLM <- function( list( model = model, formula = form, + converged = model$converged, + boundary = model$boundary, AICc = AICc ) ) @@ -527,6 +621,8 @@ trainGLM <- function( out <- data.frame( formula = form, + converged = model$converged, + boundary = model$boundary, AICc = AICc ) From d326e5f017bde2b601ba9da348f3e198428c8977 Mon Sep 17 00:00:00 2001 From: "Adam B. Smith" Date: Sat, 7 Dec 2024 18:59:00 -0600 Subject: [PATCH 2/6] Fix bug re fadctors, removes unconverged --- R/trainGlm.r | 46 ++++--- R/trainNs.r | 349 +++++++++++++++++++++++++++++++----------------- man/trainGlm.Rd | 8 +- man/trainNs.Rd | 8 +- 4 files changed, 269 insertions(+), 142 deletions(-) diff --git a/R/trainGlm.r b/R/trainGlm.r index 1ef27fd..37fb071 100644 --- a/R/trainGlm.r +++ b/R/trainGlm.r @@ -74,6 +74,7 @@ trainGLM <- function( if (FALSE) { resp <- 'presBg' + construct <- TRUE select <- TRUE quadratic <- TRUE @@ -90,11 +91,6 @@ trainGLM <- function( out <- 'model' cores <- 1 verbose <- TRUE - - # speed <- TRUE - # method <- 'eigen' - - # speed <- FALSE method <- 'glm.fit' } @@ -167,8 +163,9 @@ trainGLM <- function( ### create vector of terms terms <- preds - if (quadratic) terms <- c(terms, .makeQuadsMarginality(preds=preds, n=n, presPerTermInitial=presPerTermInitial)) - if (interaction) terms <- c(terms, .makeIAsMarginality(preds=preds, n=n, presPerTermInitial=presPerTermInitial)) + factors <- sapply(data[ , preds, drop = FALSE], is.factor) + if (quadratic) terms <- c(terms, .makeQuadsMarginality(preds = preds, n = n, presPerTermInitial = presPerTermInitial, factors = factors)) + if (interaction) terms <- c(terms, .makeIAsMarginality(preds = preds, n = n, presPerTermInitial = presPerTermInitial)) ## term-by-term model construction ################################## @@ -253,11 +250,13 @@ trainGLM <- function( } thisForm <- paste0(resp, ' ~ 1 + ', form) - numTerms <- lengths(regmatches(thisForm, gregexpr('\\+', thisForm))) + 1 - start <- rep(0, numTerms) + thisForm <- stats::as.formula(thisForm) + + mm <- stats::model.matrix(thisForm, data) + start <- rep(0, ncol(mm)) model <- suppressWarnings(stats::glm( - formula = stats::as.formula(thisForm), + formula = thisForm, family = family, data = data, method = method, @@ -479,11 +478,13 @@ trainGLM <- function( form <- unique(form) form <- paste(form, collapse = ' + ') thisForm <- paste0(resp, ' ~ 1 + ', form) - numTerms <- lengths(regmatches(thisForm, gregexpr('\\+', thisForm))) + 1 - start <- rep(0, numTerms) + thisForm <- stats::as.formula(thisForm) + + mm <- stats::model.matrix(thisForm, data) + start <- rep(0, ncol(mm)) model <- suppressWarnings(stats::glm( - formula = stats::as.formula(thisForm), + formula = thisForm, family = family, data = data, method = method, @@ -589,8 +590,9 @@ trainGLM <- function( } } thisForm <- paste0(resp, ' ~ ', form) - numTerms <- lengths(regmatches(thisForm, gregexpr('\\+', thisForm))) + 1 - start <- rep(0, numTerms) + thisForm <- stats::as.formula(thisForm) + mm <- stats::model.matrix(thisForm, data) + start <- rep(0, ncol(mm)) model <- suppressWarnings(stats::glm( formula = stats::as.formula(thisForm), @@ -632,11 +634,21 @@ trainGLM <- function( } # make vector of quadratic terms respecting marginality -.makeQuadsMarginality <- function(preds, n, presPerTermInitial) { +.makeQuadsMarginality <- function( + preds, # vector of predictor names + n, # sample size + presPerTermInitial, # number of presences per term + factors # logical: is each predictor a factor? +) { quads <- character() if (n >= 2 * presPerTermInitial) { - for (i in seq_along(preds)) quads <- c(quads, paste0(preds[i], ' + I(', preds[i], '^2)')) + for (i in seq_along(preds)) { + pred <- preds[i] + if (!factors[[pred]]) { + quads <- c(quads, paste0(pred, ' + I(', pred, '^2)')) + } + } } quads diff --git a/R/trainNs.r b/R/trainNs.r index 762e9f5..56a493f 100644 --- a/R/trainNs.r +++ b/R/trainNs.r @@ -11,7 +11,7 @@ #' @param interaction If \code{TRUE} (default), include two-way interaction terms. #' @param interceptOnly If \code{TRUE} (default) and model selection is enabled, then include an intercept-only model. #' @param presPerTermFinal Minimum number of presence sites per term in initial starting model. -#' @param maxTerms Maximum number of terms to be used in any model, not including the intercept (default is 8). Used only if \code{construct} is \code{TRUE}. +#' @param maxTerms Maximum number of terms to be used in any model, not including the intercept (default is 8). #' @param w Weights. Any of: #' \itemize{ #' \item \code{TRUE}: Causes the total weight of presences to equal the total weight of absences (if \code{family='binomial'}) @@ -20,6 +20,8 @@ #' \item The name of the column in \code{data} that contains site weights. #' } #' @param family Name of family for data error structure (see \code{\link[stats]{family}}). +#' @param removeInvalid Logical. If \code{TRUE} (default), remove models that either did not converge or have parameter estimates near the boundaries (usually negative or positive infinity). +#' @param failIfInvalid Logical. If \code{TRUE} (default), and the "full" model either does not converge or has parameters near the boundary, then the function will fail. If \code{FALSE}, then return \code{NULL} in this case. #' @param out Character vector. One or more values: #' \itemize{ #' \item \code{'model'}: Model with the lowest AICc. @@ -50,12 +52,35 @@ trainNS <- function( maxTerms = 8, w = TRUE, family = 'binomial', + removeInvalid = TRUE, + failIfInvalid = TRUE, out = 'model', cores = 1, verbose = FALSE, ... ) { + if (FALSE) { + + resp <- 'presBg' + + scale <- TRUE + df <- 1:4 + interaction <- TRUE + interceptOnly <- TRUE + method <- 'glm.fit' + presPerTermFinal <- 10 + maxTerms <- 8 + w <- TRUE + family <- 'binomial' + removeInvalid <- TRUE + failIfInvalid <- TRUE + out <- 'model' + cores <- 1 + verbose <- TRUE + + } + ########### ## setup ## ########### @@ -101,32 +126,44 @@ trainNS <- function( paths <- .libPaths() # need to pass this to avoid "object '.doSnowGlobals' not found" error!!! mcOptions <- list(preschedule = TRUE, set.seed = TRUE, silent = verbose) - ### make list of candidate model terms - ###################################### - n <- if (family %in% c('binomial', 'quasibinomial')) { sum(data[ , resp, drop=TRUE]) } else { nrow(data) } - ### create vector of terms + ### model selection + ################### + + ### make list of candidate model terms + terms <- character() + factors <- sapply(data[ , preds, drop = FALSE], is.factor) + names(factors) <- preds + componentTerms <- character() + + # univariate terms for (thisPred in preds) { - for (thisDf in df) { - - term <- if (!inherits(data[ , thisPred], 'factor')) { - paste0('splines::ns(', thisPred, ', df=', thisDf, ')') - } else { - thisPred - } + if (factors[thisPred]) { + terms <- c(terms, thisPred) + componentTerms <- c(componentTerms, thisPred) + } else { + + for (thisDf in df) { - terms <- c(terms, term) + term <- paste0('splines::ns(', thisPred, ', df=', thisDf, ')') + terms <- c(terms, term) + componentTerms <- c(componentTerms, thisPred) + + } + } } + numTerms <- rep(1, length(terms)) + # interaction terms if (interaction & length(preds) > 1L & n >= 2 * presPerTermFinal) { @@ -137,87 +174,138 @@ trainNS <- function( pred1 <- predCombos[[i]][1] pred2 <- predCombos[[i]][2] - terms <- c(terms, paste0('splines::ns(', pred1, ' * ', pred2, ', df=', df, ')')) + if (factors[[pred1]] & !factors[[pred2]]) { + newTerm <- paste0(pred1, '* splines::ns(', pred2, ', df=', df, ')') + } else if (!factors[[pred1]] & factors[[pred2]]) { + newTerm <- paste0(pred2, '* splines::ns(', pred1, ', df=', df, ')') + } else { + newTerm <- paste0('splines::ns(', pred1, ' * ', pred2, ', df=', df, ')') + } + + componentTerms <- c(componentTerms, rep(paste(pred1, pred2, collapse = ' '), length(df))) + numTerms <- c(numTerms, rep(2, length(df))) + terms <- c(terms, newTerm) } } # if more than one term - - terms <- sort(terms) - ### model selection - ################### + ### select best set of terms for full model + ########################################### - ### make table of all possible terms - - # single-predictor terms - candidates <- list(x = paste0('splines::ns(', preds[1L], ', df=', df, ')')) - candidates$x <- c(candidates$x, NA) + forms <- terms + # if (interceptOnly) forms <- c(forms, '1') - if (length(preds) > 1L) { - for (i in 2L:length(preds)) { - candidates <- c( - candidates, - list(x = c(paste0('splines::ns(', preds[i], ', df=', df, ')'), NA)) - ) + if (verbose) omnibus::say('Evaluating simple models with each candidate term...') + work <- foreach::foreach( + i = seq_along(forms), + .options.multicore = mcOptions, + .combine = 'rbind', + .inorder = FALSE, + .export = c('.trainNsWorker') + ) %makeWork% { + .trainNsWorker( + i = i, + forms = forms, + numTerms = numTerms, + componentTerms = componentTerms, + data = data, + resp = resp, + family = family, + method = method, + w = w, + insertIntercept = FALSE, + paths = paths, + modelOut = FALSE, + ... + ) + } + + if (removeInvalid) { + + bads <- which(!work$converged | work$boundary) + if (length(bads) > 0) { + work <- work[-bads, , drop = FALSE] + + if (nrow(work) == 0) { + msg <- 'No single-term models converged or all models had parameter estimates near the boundary.' + if (failIfInvalid) { + stop(msg) + } else { + warning(msg) + return(NULL) + } + } } - } - names(candidates) <- paste0('pred', seq_along(preds)) - - # add interaction terms - if (interaction & length(preds) > 1L & n >= 2 * presPerTermFinal) { - - for (countPred1 in 1L:(length(preds) - 1L)) { # for each predictor test two-variable terms - pred1 <- preds[countPred1] + work <- work[order(work$AICc), , drop = FALSE] + rownames(work) <- NULL - for (countPred2 in 2L:length(preds)) { # for each second predictor test two-variable terms + if (verbose) { + omnibus::say('Term-by-term evaluation:', pre = 1) + print(work[ , c('term', 'converged', 'boundary', 'AICc')]) + utils::flush.console() + } - pred2 <- preds[countPred2] - candidates$x <- c(paste0('splines::ns(', preds[countPred1], ' * ', preds[countPred2], ', df=', df, ')'), NA) - names(candidates)[length(candidates)] <- paste0('pred', length(candidates)) - - } # next second term - - } # next first term + ### model construction + ###################### + + topTerms <- work$term[1L] + totalTerms <- work$numTerms[1L] + + termsSoFar <- work$componentTerms[1L] - } # if more than one term + i <- 2L + while (i <= nrow(work) && (n >= (totalTerms + numTerms[i]) * presPerTermFinal & totalTerms < maxTerms + numTerms[i])) { - candidates <- expand.grid(candidates, stringsAsFactors = FALSE) - - numTerms <- rowSums(!is.na(candidates)) - candidates <- candidates[numTerms > 0, , drop=FALSE] + if (!(work$componentTerms[i] %in% termsSoFar)) { + topTerms <- c(topTerms, work$term[i]) + totalTerms <- totalTerms + work$numTerms[i] + termsSoFar <- c(termsSoFar, work$componentTerms[i]) + } - numTerms <- rowSums(!is.na(candidates)) - candidates <- candidates[numTerms <= maxTerms, , drop=FALSE] + i <- i + 1L - numTerms <- rowSums(!is.na(candidates)) - candidates <- candidates[n >= numTerms * presPerTermFinal, , drop=FALSE] - - forms <- rep(NA, nrow(candidates)) - for (i in 1L:nrow(candidates)) { - - terms <- which(!is.na(unlist(candidates[i , ]))) - terms <- unlist(candidates[i, terms]) - if (length(terms) > 1L) terms <- paste(terms, collapse = ' + ') - forms[i] <- terms - } + + formsNA <- list() + for (i in seq_along(topTerms)) { + formsNA[[i]] <- c(topTerms[i], NA_character_) + } + + formGrid <- expand.grid(formsNA, stringsAsFactors = FALSE) + + forms <- character() + for (i in 1:nrow(formGrid)) { - forms <- paste0('1 + ', forms) - if (interceptOnly) forms <- c(forms, '1') + terms <- unlist(formGrid[i, , drop = TRUE]) + terms <- terms[!is.na(terms)] + if (length(terms) > 0) { + form <- paste(terms, collapse = ' + ') + forms[i] <- form + } + } + + + if (verbose) omnibus::say('Evaluating all possible models constructed from full model...', pre = 1) + + wantModels <- any(c('models', 'model') %in% out) + combine <- if (wantModels) { 'c' } else { 'rbind' } + work <- foreach::foreach( i = seq_along(forms), .options.multicore = mcOptions, - .combine = 'c', + .combine = combine, .inorder = FALSE, .export = c('.trainNsWorker') ) %makeWork% { .trainNsWorker( i = i, forms = forms, + numTerms = numTerms, + componentTerms = componentTerms, data = data, resp = resp, family = family, @@ -225,59 +313,65 @@ trainNS <- function( w = w, insertIntercept = FALSE, paths = paths, - modelOut = ('models' %in% out | 'model' %in% out), + modelOut = wantModels, ... ) } - - # if (cores > 1L) parallel::stopCluster(cl) - # tuning table - tuning <- data.frame( - model = work[[1L]]$formula, - AICc = work[[1L]]$AICc - ) + if (wantModels) models <- list() + tuning <- data.frame() + for (i in seq_along(work)) { - if (length(work) > 1L) { - for (i in 2L:length(work)) { - - tuning <- rbind( - tuning, - data.frame( - model = work[[i]]$formula, - AICc = work[[i]]$AICc - ) + if (wantModels) { + models[[i]] <- work[[i]]$model + models[[i]]$scales <- scales + } + + tuning <- rbind( + tuning, + data.frame( + model = work[[i]]$formula, + converged = work[[i]]$model$converged, + boundary = work[[i]]$model$boundary, + AICc = work[[i]]$AICc ) + ) + + } + + if (removeInvalid) { + + bads <- which(!tuning$converged | tuning$boundary) + if (length(bads) > 0) { + tuning <- work[-bads, , drop = FALSE] + + if (nrow(work) == 0) { + msg <- 'No single-term models converged or all models had parameter estimates near the boundary.' + if (failIfInvalid) { + stop(msg) + } else { + warning(msg) + return(NULL) + } + } + + if (wantModels) models <- models[-bads] + } + } - - if (!is.na(scale)) { - if (scale) for (i in seq_along(work)) work[[i]]$model$scale <- scales - } - + bestOrder <- order(tuning$AICc) - if ('model' %in% out) model <- work[[bestOrder[1L]]]$model - if ('models' %in% out) { - models <- list() - models[[1]] <- work[[1L]]$model - for (i in 2L:length(work)) models[[i]] <- work[[i]]$model - models <- models[bestOrder] - } tuning <- tuning[bestOrder, , drop = FALSE] rownames(tuning) <- NULL - - if (verbose) { - - omnibus::say('Model selection:', level=2) - print(tuning) - utils::flush.console() - + if (wantModels) { + models <- models[bestOrder] + model <- models[[1L]] } ### return ########## - if (length(out) > 1L) { output <- list() if ('models' %in% out) output$models <- models @@ -294,13 +388,14 @@ trainNS <- function( } -################# -### train GAM ### -################# +### train NS +############ .trainNsWorker <- function( i, - forms, # formulae (without LHS and maybe without intercept) + forms, # formulae (without LHS and intercept) + numTerms, # vector of number of predictors in each term + componentTerms, # character vector of predictors in each term data, resp, family, @@ -315,24 +410,25 @@ trainNS <- function( # need to call this to avoid "object '.doSnowGlobals' not found" error!!! .libPaths(paths) - form <- forms[i] - if (insertIntercept) { - form <- if (form == '') { - '1' - } else { - paste('1', form, sep=' + ') - } - } - thisForm <- paste0(resp, ' ~ ', form) - - model <- stats::glm( - formula = stats::as.formula(thisForm), + # form <- forms[i] + # if (insertIntercept) { + # form <- if (form == '') { + # '1' + # } else { + # paste('1', form, sep=' + ') + # } + # } + thisForm <- paste0(resp, ' ~ ', forms[i]) + thisForm <- stats::as.formula(thisForm) + + model <- suppressWarnings(stats::glm( + formula = thisForm, data = data, family = family, method = method, weights = w, ... - ) + )) AICc <- AICcmodavg::AICc(model) @@ -342,7 +438,11 @@ trainNS <- function( list( list( model = model, - formula = form, + formula = forms[i], + numTerms = numTerms[i], + componentTerms = componentTerms[i], + converged = model$converged, + boundary = model$boundary, AICc = AICc ) ) @@ -350,12 +450,15 @@ trainNS <- function( } else { data.frame( - formula = form, + term = forms[i], + componentTerms = componentTerms[i], + numTerms = numTerms[i], + converged = model$converged, + boundary = model$boundary, AICc = AICc ) } - out } diff --git a/man/trainGlm.Rd b/man/trainGlm.Rd index ae0c887..f693188 100644 --- a/man/trainGlm.Rd +++ b/man/trainGlm.Rd @@ -20,6 +20,8 @@ trainGLM( maxTerms = 8, w = TRUE, family = stats::binomial(), + removeInvalid = TRUE, + failIfInvalid = TRUE, out = "model", cores = 1, verbose = FALSE, @@ -63,6 +65,10 @@ trainGLM( \item{family}{Name of family for data error structure (see \code{\link[stats]{family}}). Default is to use the 'binomial' family.} +\item{removeInvalid}{Logical. If \code{TRUE} (default), remove models that either did not converge or have parameter estimates near the boundaries (usually negative or positive infinity). If you run this function with `construct = TRUE` (i.e., construct a "full" model from the best "small" models), then any small model that either did not converge or had parameters that are near the boundary (usually negative or positive infinity) are removed from consideration as terms in "full" model.} + +\item{failIfInvalid}{Logical. If \code{TRUE} (default), and the "full" model either does not converge or has parameters near the boundary, then the function will fail. If \code{FALSE}, then return \code{NULL} in this case.} + \item{out}{Character vector. One or more values: \itemize{ \item \code{'model'}: Model with the lowest AICc. @@ -77,7 +83,7 @@ trainGLM( \item{...}{Arguments to pass to \code{glm}.} } \value{ -The object that is returned depends on the value of the \code{out} argument. It can be a model object, a data frame, a list of models, or a list of all two or more of these. If \code{scale} is \code{TRUE}, any model object will also have an element named \code{$scale}, which contains the means and standard deviations for predictors that are not factors. +The object that is returned depends on the value of the \code{out} argument. It can be a model object, a data frame, a list of models, or a list of all two or more of these. If \code{scale} is \code{TRUE}, any model object will also have an element named \code{$scale}, which contains the means and standard deviations for predictors that are not factors. The data frame reports the AICc for all of the models, sorted by best to worst. The \code{converged} column indicates whether the model converged ("\code{TRUE}" is good), and the \code{boundary} column whether the model parameters are near the boundary (usually, negative or positive infinity; "\code{FALSE}" is good). } \description{ This function constructs a generalized linear model. By default, the model is constructed in a two-stage process. First, the "construct" phase generates a series of simple models with univariate, quadratic, or 2-way-interaction terms. These simple models are then ranked based on their AICc. Second, the "select" phase creates a "full" model from the simple models such that there is at least \code{presPerTermInitial} presences (if the response is binary) or data rows (if not) for each coefficient to be estimated (not counting the intercept). Finally, it selects the best model using AICc from all possible subsets of this "full" model, while respecting marginality (i.e., all lower-order terms of higher-order terms appear in the model). diff --git a/man/trainNs.Rd b/man/trainNs.Rd index f40dc61..ee981ba 100644 --- a/man/trainNs.Rd +++ b/man/trainNs.Rd @@ -17,6 +17,8 @@ trainNS( maxTerms = 8, w = TRUE, family = "binomial", + removeInvalid = TRUE, + failIfInvalid = TRUE, out = "model", cores = 1, verbose = FALSE, @@ -42,7 +44,7 @@ trainNS( \item{presPerTermFinal}{Minimum number of presence sites per term in initial starting model.} -\item{maxTerms}{Maximum number of terms to be used in any model, not including the intercept (default is 8). Used only if \code{construct} is \code{TRUE}.} +\item{maxTerms}{Maximum number of terms to be used in any model, not including the intercept (default is 8).} \item{w}{Weights. Any of: \itemize{ @@ -54,6 +56,10 @@ trainNS( \item{family}{Name of family for data error structure (see \code{\link[stats]{family}}).} +\item{removeInvalid}{Logical. If \code{TRUE} (default), remove models that either did not converge or have parameter estimates near the boundaries (usually negative or positive infinity).} + +\item{failIfInvalid}{Logical. If \code{TRUE} (default), and the "full" model either does not converge or has parameters near the boundary, then the function will fail. If \code{FALSE}, then return \code{NULL} in this case.} + \item{out}{Character vector. One or more values: \itemize{ \item \code{'model'}: Model with the lowest AICc. From 153983368c4724b15f32e800f0e9bfd15105a431 Mon Sep 17 00:00:00 2001 From: "Adam B. Smith" Date: Sat, 7 Dec 2024 18:59:04 -0600 Subject: [PATCH 3/6] Update enmSdmX_workspace.code-workspace --- enmSdmX_workspace.code-workspace | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/enmSdmX_workspace.code-workspace b/enmSdmX_workspace.code-workspace index f151917..0e314d8 100644 --- a/enmSdmX_workspace.code-workspace +++ b/enmSdmX_workspace.code-workspace @@ -24,6 +24,7 @@ "titleBar.inactiveBackground": "#1857a499", "titleBar.inactiveForeground": "#e7e7e799" }, - "peacock.color": "#1857a4" + "peacock.color": "#1857a4", + "git.ignoreLimitWarning": true } } \ No newline at end of file From 617537dd28d8e8a035487c3dc6e11f67c75a1d58 Mon Sep 17 00:00:00 2001 From: "Adam B. Smith" Date: Sat, 7 Dec 2024 18:59:08 -0600 Subject: [PATCH 4/6] Update NEWS.md --- NEWS.md | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/NEWS.md b/NEWS.md index f556bd9..7c336b4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,9 @@ -# enmSdmX 1.1.10 2024-12-03 -- `trainNS()` works when the number of predictors is >2 and interactions are allowed between variables (thank you again, Pascal!). +# enmSdmX 1.2.10 2024-12-07 +- `trainGLM()` now allows users to remove models that did not converge or had boundary problems. The function also works for factor predictors with > 2 levels. Thank you, A.L.! +- `trainNS()` works when the number of predictors is >2 and interactions are allowed between variables. Thank you, P.T.! It also handles factors with >2 levels and allows users to remove models that did not converge or had boundary problems. Thank you, A.L.! # enmSdmX 1.1.9 2024-11-01 -- `trainGAM()` does not include interaction terms if `interaction` is `NULL`, and does not fail when interactions are included and the number of predictors is >2 (thank you, Pascal!). +- `trainGAM()` does not include interaction terms if `interaction` is `NULL`, and does not fail when interactions are included and the number of predictors is >2 (thank you, P.T.!). # enmSdmX 1.1.8 2024-10-02 - `modelSize()` can now tell size of a `ranger` random forest. From d5902c411617b8fb76b25030df3fcfcb045d8eb8 Mon Sep 17 00:00:00 2001 From: "Adam B. Smith" Date: Sat, 7 Dec 2024 18:59:12 -0600 Subject: [PATCH 5/6] Update aaa.r --- R/aaa.r | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/R/aaa.r b/R/aaa.r index 189981c..5e9309b 100644 --- a/R/aaa.r +++ b/R/aaa.r @@ -6,8 +6,12 @@ ver <- read.dcf(file=system.file('DESCRIPTION', package=pkg), fields='Version') packageStartupMessage(paste0('This is ', pkg, ' ', ver, '.')) - packageStartupMessage(paste0('Back-incompatible changes starting starting with enmSdmX version 1.1.1:')) + + packageStartupMessage(paste0('Back-incompatible changes starting starting with enmSdmX version 1.1.1 (2023-06-11):')) packageStartupMessage(paste0('* trainRF() replaces use of the randomForest package with the faster ranger package, which produces random forests that are statistically equivalent.')) packageStartupMessage(paste0('* geoFold() uses ', dQuote('complete'), ' clustering by default (was ', dQuote('single'), ').')) + packageStartupMessage(paste0('Back-incompatible changes starting starting with enmSdmX version 1.2.10 (2024-12-07):')) + packageStartupMessage(paste0('* By default, trainGLM() and trainNS() remove unconverged models and models with boundary issues.')) + } From b65c034fda7f61e35eece7624e265cb36ff60056 Mon Sep 17 00:00:00 2001 From: "Adam B. Smith" Date: Sat, 7 Dec 2024 18:59:15 -0600 Subject: [PATCH 6/6] Update DESCRIPTION --- DESCRIPTION | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index dfe3ecf..1897634 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: enmSdmX Type: Package Title: Species Distribution Modeling and Ecological Niche Modeling -Version: 1.1.10 -Date: 2024-12-03 +Version: 1.2.10 +Date: 2024-12-07 Authors@R: c( person(