manual merge develop --> main

.Rbuildignore

@@ -33,3 +33,4 @@ deprecated
 .github/
 ^\.github$
 ^CRAN-RELEASE$
+.covrignore

.covrignore

@@ -0,0 +1,10 @@
+src/cyclops/io/*InputReader.*
+src/cyclops/drivers/Bootstrap*
+src/cyclops/drivers/Hierarchy*
+src/cyclops/drivers/GridSearch*
+src/cyclops/drivers/Proportion*
+src/cyclops/drivers/CrossValidationSelector.h
+src/cyclops/drivers/AbstractSelector.h
+R/cyclops.R
+R/PackageMaintenance.R
+

.gitignore

@@ -30,3 +30,5 @@ extras
 CRAN-RELEASE
 .Rdata
 .DS_Store
+*.gcno
+*.gcda

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: Cyclops
 Type: Package
 Title: Cyclic Coordinate Descent for Logistic, Poisson and Survival Analysis
-Version: 3.1.2
+Version: 3.1.2.9999
 Authors@R: c(
     person("Marc A.", "Suchard", email = "[email protected]", role = c("aut","cre")),
     person("Martijn J.", "Schuemie", role = "aut"),
@@ -44,6 +44,7 @@ LinkingTo: Rcpp,
     RcppEigen (>= 0.3.2)
 Suggests:
     testthat,
+    readr,
     MASS,
     gnm,
     ggplot2,

NEWS.md

@@ -1,3 +1,10 @@
+Cyclops v3.1.3
+
+Changes:
+
+1. fixed likelihood profiling when non-convex due to numerical instability
+2. fixed parsing of 64-bit covariate IDs
+
 Cyclops v3.1.2
 ==============
 

R/ModelFit.R

@@ -32,6 +32,7 @@
 #' @param returnEstimates Logical, return regression coefficient estimates in Cyclops model fit object
 #' @param startingCoefficients Vector of starting values for optimization
 #' @param fixedCoefficients Vector of booleans indicating if coefficient should be fix
+#' @param warnings Logical, report regularization warnings
 #' @param computeDevice String: Name of compute device to employ; defaults to \code{"native"} C++ on CPU
 #'
 #' @return
@@ -70,6 +71,7 @@ fitCyclopsModel <- function(cyclopsData,
                             returnEstimates = TRUE,
                             startingCoefficients = NULL,
                             fixedCoefficients = NULL,
+                            warnings = TRUE,
 							computeDevice = "native") {
 
     # Delegate to control$setHook if exists
@@ -104,7 +106,7 @@ fitCyclopsModel <- function(cyclopsData,
 
     if (!is.null(prior$setHook)) {
 
-        prior$setHook(cyclopsData) # Call-back
+        prior$setHook(cyclopsData, warnings) # Call-back
 
     } else {
         prior$exclude <- .checkCovariates(cyclopsData, prior$exclude)
@@ -132,7 +134,7 @@ fitCyclopsModel <- function(cyclopsData,
                     warn <- TRUE
                 }
             }
-            if (warn) {
+            if (warn && warnings) {
                 warning("Excluding intercept from regularization")
             }
         }
@@ -167,6 +169,22 @@ fitCyclopsModel <- function(cyclopsData,
             }
         }
 
+        if (any(prior$priorType == "jeffreys")) {
+            if (Cyclops::getNumberOfCovariates(cyclopsData) > 1) {
+                stop("Jeffreys prior is currently only implemented for 1 covariate")
+            }
+
+            covariate <- Cyclops::getCovariateIds(cyclopsData)
+            if (Cyclops::getCovariateTypes(cyclopsData, covariate) != "indicator") {
+                count <- reduce(cyclopsData, covariate, power = 0)
+                sum <- reduce(cyclopsData, covariate, power = 1)
+                mean <- sum / count
+                if (!(mean == 0.0 || mean == 1.0)) {
+                    stop("Jeffreys prior is currently only implemented for indicator covariates")
+                }
+            }
+        }
+
         .cyclopsSetPrior(cyclopsData$cyclopsInterfacePtr, prior$priorType, prior$variance,
                          prior$exclude, graph, neighborhood)
     }
@@ -247,7 +265,7 @@ fitCyclopsModel <- function(cyclopsData,
     }
 
     if (!is.null(cyclopsData$censorWeights)) {
-        if (cyclopsData$modelType != 'fgr') {
+        if (cyclopsData$modelType != 'fgr' && warnings) {
             warning(paste0("modelType = '", cyclopsData$modelType, "' does not use censorWeights. These weights will not be passed further."))
         }
         if (length(cyclopsData$censorWeights) != getNumberOfRows(cyclopsData)) {
@@ -273,6 +291,24 @@ fitCyclopsModel <- function(cyclopsData,
         fit <- .cyclopsFitModel(cyclopsData$cyclopsInterfacePtr)
     }
 
+    if (fit$return_flag == "POOR_BLR_STEP" && control$convergenceType == "gradient") {
+
+        if (warnings) {
+            warning("BLR convergence criterion failed; coefficient may be infinite")
+        }
+
+        control$convergenceType <- "lange"
+        return(fitCyclopsModel(cyclopsData = cyclopsData,
+                               prior = prior,
+                               control = control,
+                               weights = weights,
+                               forceNewObject = forceNewObject,
+                               returnEstimates = returnEstimates,
+                               startingCoefficients = startingCoefficients,
+                               fixedCoefficients = fixedCoefficients,
+                               computeDevice = computeDevice))
+    }
+
     if (returnEstimates) {
         estimates <- .cyclopsLogModel(cyclopsData$cyclopsInterfacePtr)
         fit <- c(fit, estimates)
@@ -295,12 +331,14 @@ fitCyclopsModel <- function(cyclopsData,
 .checkCovariates <- function(cyclopsData, covariates) {
     if (!is.null(covariates)) {
         saved <- covariates
+
+        indices <- NULL
+
         if (inherits(covariates, "character")) {
             # Try to match names
             indices <- match(covariates, cyclopsData$coefficientNames)
             covariates <- getCovariateIds(cyclopsData)[indices]
         }
-        # covariates = as.numeric(covariates)
 
         if (!bit64::is.integer64(covariates)) {
             covariates <- bit64::as.integer64(covariates)
@@ -309,6 +347,8 @@ fitCyclopsModel <- function(cyclopsData,
         if (any(is.na(covariates))) {
             stop("Unable to match all covariates: ", paste(saved, collapse = ", "))
         }
+
+        attr(covariates, "indices") <- indices
     }
     covariates
 }
@@ -582,7 +622,7 @@ createPrior <- function(priorType,
                         neighborhood = NULL,
                         useCrossValidation = FALSE,
                         forceIntercept = FALSE) {
-    validNames = c("none", "laplace","normal", "barupdate", "hierarchical")
+    validNames = c("none", "laplace","normal", "barupdate", "hierarchical", "jeffreys")
     stopifnot(priorType %in% validNames)
     if (!is.null(exclude)) {
         if (!inherits(exclude, "character") &&
@@ -660,7 +700,7 @@ getCyclopsPredictiveLogLikelihood <- function(object, weights) {
     }
     # TODO Remove code duplication with weights section of fitCyclopsModel
 
-    .cyclopsGetPredictiveLogLikelihood(object$cyclopsData$cyclopsInterfacePtr, weights)
+    .cyclopsGetNewPredictiveLogLikelihood(object$cyclopsData$cyclopsInterfacePtr, weights)
 }
 
 #' @title Get cross-validation information from a Cyclops model fit
@@ -690,7 +730,7 @@ getCrossValidationInfo <- function(object) {
             control$seed <- as.integer(Sys.time())
         }
 
-        if (is.null(control$algoritm) || is.na(control$algorithm)) { # Provide backwards compatibility
+        if (is.null(control$algorithm) || is.na(control$algorithm)) { # Provide backwards compatibility
             control$algorithm <- "ccd"
         }
 
@@ -770,6 +810,7 @@ confint.cyclopsFit <- function(object, parm, level = 0.95, #control,
                                rescale = FALSE, ...) {
     .checkInterface(object$cyclopsData, testOnly = TRUE)
     #.setControl(object$cyclopsData$cyclopsInterfacePtr, control)
+
     parm <- .checkCovariates(object$cyclopsData, parm)
     if (level < 0.01 || level > 0.99) {
         stop("level must be between 0 and 1")
@@ -787,12 +828,15 @@ confint.cyclopsFit <- function(object, parm, level = 0.95, #control,
                                  threads, threshold,
                                  overrideNoRegularization,
                                  includePenalty)
+
+    indices <- match(parm, getCovariateIds(object$cyclopsData))
+
     if (!is.null(object$scale) && rescale) {
-        prof$lower <- prof$lower * object$scale[as.integer(parm)]
-        prof$upper <- prof$upper * object$scale[as.integer(parm)]
+        prof$lower <- prof$lower * object$scale[indices]
+        prof$upper <- prof$upper * object$scale[indices]
     }
     prof <- as.matrix(as.data.frame(prof))
-    rownames(prof) <- object$coefficientNames[as.integer(parm)]
+    rownames(prof) <- object$coefficientNames[indices]
     qs <- c((1 - level) / 2, 1 - (1 - level) / 2) * 100
     colnames(prof)[2:3] <- paste(sprintf("%.1f", qs), "%")
 

R/NewDataConversion.R

@@ -309,13 +309,12 @@ convertToCyclopsData.tbl_dbi <- function(outcomes,
             (select(outcomes, .data$y) %>% distinct() %>% count() %>% collect() > 2)) {
             stop("Cox model only accepts one outcome type")
         }
-
-        outcomes <- outcomes %>%
-            arrange(.data$stratumId, desc(.data$time), .data$y, .data$rowId)
         if (!"time" %in% colnames(covariates)) {
             covariates <- covariates %>%
                 inner_join(select(outcomes, .data$rowId, .data$time, .data$y), by = "rowId")
         }
+        outcomes <- outcomes %>%
+            arrange(.data$stratumId, desc(.data$time), .data$y, .data$rowId)
         covariates <- covariates %>%
             arrange(.data$covariateId, .data$stratumId, desc(.data$time), .data$y, .data$rowId)
     }

R/ParameterizedPrior.R

@@ -41,7 +41,7 @@ createParameterizedPrior <- function(priorType,
         stop("Cannot perform cross validation with a flat prior")
     }
 
-    setHook <- function(cyclopsData) {
+    setHook <- function(cyclopsData, warnings) {
         # closure to capture arguments
         if (length(priorType) > 1) {
             if (length(priorType) != getNumberOfCovariates(cyclopsData)) {
@@ -51,7 +51,9 @@ createParameterizedPrior <- function(priorType,
 
         if (priorType[1] != "none" && .cyclopsGetHasIntercept(cyclopsData) && !forceIntercept) {
             priorType[1] <- "none"
-            warning("Excluding intercept from regularization")
+            if (warnings) {
+                warning("Excluding intercept from regularization")
+            }
         }
 
         .cyclopsSetParameterizedPrior(cyclopsData$cyclopsInterfacePtr,

R/RcppExports.R

@@ -41,10 +41,6 @@
     invisible(.Call(`_Cyclops_cyclopsSetCensorWeights`, inRcppCcdInterface, weights))
 }
 
-.cyclopsGetPredictiveLogLikelihood <- function(inRcppCcdInterface, weights) {
-    .Call(`_Cyclops_cyclopsGetPredictiveLogLikelihood`, inRcppCcdInterface, weights)
-}
-
 .cyclopsGetNewPredictiveLogLikelihood <- function(inRcppCcdInterface, weights) {
     .Call(`_Cyclops_cyclopsGetNewPredictiveLogLikelihood`, inRcppCcdInterface, weights)
 }
@@ -53,6 +49,10 @@
     .Call(`_Cyclops_cyclopsGetLogLikelihood`, inRcppCcdInterface)
 }
 
+.cyclopsLogResults <- function(inRcppCcdInterface, fileName, withASE) {
+    invisible(.Call(`_Cyclops_cyclopsLogResult`, inRcppCcdInterface, fileName, withASE))
+}
+
 .cyclopsGetFisherInformation <- function(inRcppCcdInterface, sexpBitCovariates) {
     .Call(`_Cyclops_cyclopsGetFisherInformation`, inRcppCcdInterface, sexpBitCovariates)
 }
@@ -109,18 +109,6 @@
     .Call(`_Cyclops_isSortedVectorList`, vectorList, ascending)
 }
 
-#' @title Print row identifiers
-#'
-#' @description
-#' \code{printCcdRowIds} return the row identifiers in a Cyclops data object
-#'
-#' @param object    A Cyclops data object
-#'
-#' @keywords internal
-printCyclopsRowIds <- function(object) {
-    invisible(.Call(`_Cyclops_cyclopsPrintRowIds`, object))
-}
-
 .isRcppPtrNull <- function(x) {
     .Call(`_Cyclops_isRcppPtrNull`, x)
 }