change confint behavior when not using convergence

OHDSI · Feb 18, 2025 · 740a8c1 · 740a8c1
1 parent 699e5c3
commit 740a8c1
Show file tree

Hide file tree

Showing 7 changed files with 74 additions and 8 deletions.
diff --git a/Cyclops.Rproj b/Cyclops.Rproj
@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: 5004ae26-a9a8-4fe8-8392-938892bbd248
 
 RestoreWorkspace: No
 SaveWorkspace: No

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: Cyclops
 Type: Package
 Title: Cyclic Coordinate Descent for Logistic, Poisson and Survival Analysis
-Version: 3.5.0
+Version: 3.5.1
 Authors@R: c(
     person("Marc A.", "Suchard", email = "[email protected]", role = c("aut","cre")),
     person("Martijn J.", "Schuemie", role = "aut"),

diff --git a/NEWS.md b/NEWS.md
@@ -3,6 +3,7 @@ develop
 
 1. add Schoenfeld residual output for Cox models
 2. check for negative curvature before computing CIs
+   a. change to "lange"-convergence when needed
 3. fix `vcov` when model has an offset
 
 Cyclops v3.5.0

diff --git a/R/ModelFit.R b/R/ModelFit.R
@@ -861,14 +861,27 @@ confint.cyclopsFit <- function(object, parm, level = 0.95, #control,
 
     hessianDiagonal <- .cyclopsGetLogLikelihoodHessianDiagonal(object$cyclopsData$cyclopsInterfacePtr,
                                                                parm)
-    if (any(hessianDiagonal > maximumCurvature)) {
-        warning("Cannot estimate confidence interval for a monotonic log-likelihood function")
 
-        prof <- data.frame(covariate = savedParm,
-                           lower = rep(NA, length(parm)),
-                           upper = rep(NA, length(parm)),
-                           evaluations = rep(NA, length(parm)))
-    } else {
+    convergenceType <- .cyclopsGetConvergenceType(object$cyclopsData$cyclopsInterfacePtr)
+    prof <- NULL
+    if (convergenceType != 1 && any(hessianDiagonal >= maximumCurvature)) {
+
+        .cyclopsSetConvergenceType(object$interface, "lange")
+        newFit <- .cyclopsFitModel(cyclopsData$cyclopsInterfacePtr)
+        newHessianDiagonal <- .cyclopsGetLogLikelihoodHessianDiagonal(object$cyclopsData$cyclopsInterfacePtr,
+                                                                      parm)
+
+        if (any(newHessianDiagonal > 0)) {
+            warning("Cannot estimate confidence interval for positive-curvature log-likelihood function")
+
+            prof <- data.frame(covariate = savedParm,
+                               lower = rep(NA, length(parm)),
+                               upper = rep(NA, length(parm)),
+                               evaluations = rep(NA, length(parm)))
+        }
+    }
+
+    if (is.null(prof)) {
 
         prof <- .cyclopsProfileModel(object$cyclopsData$cyclopsInterfacePtr, parm,
                                      threads, threshold,

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -101,6 +101,14 @@
     .Call(`_Cyclops_cyclopsPredictModel`, inRcppCcdInterface)
 }
 
+.cyclopsSetConvergenceType <- function(inRcppCcdInterface, convergenceType) {
+    invisible(.Call(`_Cyclops_cyclopsSetConvergenceType`, inRcppCcdInterface, convergenceType))
+}
+
+.cyclopsGetConvergenceType <- function(inRcppCcdInterface) {
+    .Call(`_Cyclops_cyclopsGetConvergenceType`, inRcppCcdInterface)
+}
+
 .cyclopsSetControl <- function(inRcppCcdInterface, maxIterations, tolerance, convergenceType, useAutoSearch, fold, foldToCompute, lowerLimit, upperLimit, gridSteps, noiseLevel, threads, seed, resetCoefficients, startingVariance, useKKTSwindle, swindleMultipler, selectorType, initialBound, maxBoundCount, algorithm, doItAll, syncCV) {
     invisible(.Call(`_Cyclops_cyclopsSetControl`, inRcppCcdInterface, maxIterations, tolerance, convergenceType, useAutoSearch, fold, foldToCompute, lowerLimit, upperLimit, gridSteps, noiseLevel, threads, seed, resetCoefficients, startingVariance, useKKTSwindle, swindleMultipler, selectorType, initialBound, maxBoundCount, algorithm, doItAll, syncCV))
 }

diff --git a/src/RcppCyclopsInterface.cpp b/src/RcppCyclopsInterface.cpp
@@ -540,6 +540,25 @@ List cyclopsPredictModel(SEXP inRcppCcdInterface) {
 	return list;
 }
 
+// [[Rcpp::export(".cyclopsSetConvergenceType")]]
+void cyclopsSetConvergenceType(SEXP inRcppCcdInterface, const std::string& convergenceType) {
+    using namespace bsccs;
+    XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+    // Convergence control
+    CCDArguments& args = interface->getArguments();
+    args.modeFinding.convergenceType = RcppCcdInterface::parseConvergenceType(convergenceType);
+}
+
+
+// [[Rcpp::export(".cyclopsGetConvergenceType")]]
+int cyclopsGetConvergenceType(SEXP inRcppCcdInterface) {
+    using namespace bsccs;
+    XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+    // Convergence control
+    CCDArguments& args = interface->getArguments();
+    return args.modeFinding.convergenceType;
+}
+
 
 // [[Rcpp::export(".cyclopsSetControl")]]
 void cyclopsSetControl(SEXP inRcppCcdInterface,

diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -311,6 +311,28 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cyclopsSetConvergenceType
+void cyclopsSetConvergenceType(SEXP inRcppCcdInterface, const std::string& convergenceType);
+RcppExport SEXP _Cyclops_cyclopsSetConvergenceType(SEXP inRcppCcdInterfaceSEXP, SEXP convergenceTypeSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP);
+    Rcpp::traits::input_parameter< const std::string& >::type convergenceType(convergenceTypeSEXP);
+    cyclopsSetConvergenceType(inRcppCcdInterface, convergenceType);
+    return R_NilValue;
+END_RCPP
+}
+// cyclopsGetConvergenceType
+int cyclopsGetConvergenceType(SEXP inRcppCcdInterface);
+RcppExport SEXP _Cyclops_cyclopsGetConvergenceType(SEXP inRcppCcdInterfaceSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP);
+    rcpp_result_gen = Rcpp::wrap(cyclopsGetConvergenceType(inRcppCcdInterface));
+    return rcpp_result_gen;
+END_RCPP
+}
 // cyclopsSetControl
 void cyclopsSetControl(SEXP inRcppCcdInterface, int maxIterations, double tolerance, const std::string& convergenceType, bool useAutoSearch, int fold, int foldToCompute, double lowerLimit, double upperLimit, int gridSteps, const std::string& noiseLevel, int threads, int seed, bool resetCoefficients, double startingVariance, bool useKKTSwindle, int swindleMultipler, const std::string& selectorType, double initialBound, int maxBoundCount, const std::string& algorithm, bool doItAll, bool syncCV);
 RcppExport SEXP _Cyclops_cyclopsSetControl(SEXP inRcppCcdInterfaceSEXP, SEXP maxIterationsSEXP, SEXP toleranceSEXP, SEXP convergenceTypeSEXP, SEXP useAutoSearchSEXP, SEXP foldSEXP, SEXP foldToComputeSEXP, SEXP lowerLimitSEXP, SEXP upperLimitSEXP, SEXP gridStepsSEXP, SEXP noiseLevelSEXP, SEXP threadsSEXP, SEXP seedSEXP, SEXP resetCoefficientsSEXP, SEXP startingVarianceSEXP, SEXP useKKTSwindleSEXP, SEXP swindleMultiplerSEXP, SEXP selectorTypeSEXP, SEXP initialBoundSEXP, SEXP maxBoundCountSEXP, SEXP algorithmSEXP, SEXP doItAllSEXP, SEXP syncCVSEXP) {
@@ -908,6 +930,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_Cyclops_cyclopsGetProfileLikelihood", (DL_FUNC) &_Cyclops_cyclopsGetProfileLikelihood, 5},
     {"_Cyclops_cyclopsProfileModel", (DL_FUNC) &_Cyclops_cyclopsProfileModel, 6},
     {"_Cyclops_cyclopsPredictModel", (DL_FUNC) &_Cyclops_cyclopsPredictModel, 1},
+    {"_Cyclops_cyclopsSetConvergenceType", (DL_FUNC) &_Cyclops_cyclopsSetConvergenceType, 2},
+    {"_Cyclops_cyclopsGetConvergenceType", (DL_FUNC) &_Cyclops_cyclopsGetConvergenceType, 1},
     {"_Cyclops_cyclopsSetControl", (DL_FUNC) &_Cyclops_cyclopsSetControl, 23},
     {"_Cyclops_cyclopsRunCrossValidationl", (DL_FUNC) &_Cyclops_cyclopsRunCrossValidationl, 1},
     {"_Cyclops_cyclopsFitModel", (DL_FUNC) &_Cyclops_cyclopsFitModel, 1},