check for negative curvature before computing CIs

OHDSI · Jan 9, 2025 · a49a520 · a49a520
1 parent 353a0bf
commit a49a520
Show file tree

Hide file tree

Showing 9 changed files with 85 additions and 5 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -2,6 +2,7 @@ develop
 ==============
 
 1. add Schoenfeld residual output for Cox models
+2. check for negative curvature before computing CIs
 
 Cyclops v3.5.0
 ==============

diff --git a/R/ModelFit.R b/R/ModelFit.R
@@ -825,6 +825,7 @@ runBootstrap <- function(object, outFileName, treatmentId, replicates) {
 #' @param overrideNoRegularization   Logical: Enable confidence interval estimation for regularized parameters
 #' @param includePenalty    Logical: Include regularized covariate penalty in profile
 #' @param rescale   Boolean: rescale coefficients for unnormalized covariate values
+#' @param maximumCurvature Numeric: maximum curvature of log-likelihood profile required to compute CI
 #' @param ... Additional argument(s) for methods
 #'
 #' @return
@@ -838,7 +839,9 @@ runBootstrap <- function(object, outFileName, treatmentId, replicates) {
 confint.cyclopsFit <- function(object, parm, level = 0.95, #control,
                                overrideNoRegularization = FALSE,
                                includePenalty = TRUE,
-                               rescale = FALSE, ...) {
+                               rescale = FALSE,
+                               maximumCurvature = -1E-3,
+                               ...) {
     .checkInterface(object$cyclopsData, testOnly = TRUE)
     #.setControl(object$cyclopsData$cyclopsInterfacePtr, control)
 
@@ -855,6 +858,12 @@ confint.cyclopsFit <- function(object, parm, level = 0.95, #control,
         }
     }
 
+    hessianDiagonal <- .cyclopsGetLogLikelihoodHessianDiagonal(object$cyclopsData$cyclopsInterfacePtr,
+                                                               parm)
+    if (any(hessianDiagonal > maximumCurvature)) {
+        stop("Cannot estimate confidence interval for a monotonic log-likelihood function")
+    }
+
     prof <- .cyclopsProfileModel(object$cyclopsData$cyclopsInterfacePtr, parm,
                                  threads, threshold,
                                  overrideNoRegularization,

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -125,6 +125,10 @@
     .Call(`_Cyclops_cyclopsInitializeModel`, inModelData, modelType, computeDevice, computeMLE)
 }
 
+.cyclopsGetLogLikelihoodHessianDiagonal <- function(inRcppCcdInterface, sexpCovariates) {
+    .Call(`_Cyclops_cyclopsGetLogLikelihoodHessianDiagonal`, inRcppCcdInterface, sexpCovariates)
+}
+
 #' @title List available GPU devices
 #'
 #' @description

diff --git a/man/confint.cyclopsFit.Rd b/man/confint.cyclopsFit.Rd
diff --git a/man/cyclops.Rd b/man/cyclops.Rd
diff --git a/src/RcppCyclopsInterface.cpp b/src/RcppCyclopsInterface.cpp
@@ -744,6 +744,55 @@ List cyclopsInitializeModel(SEXP inModelData, const std::string& modelType, cons
 	return list;
 }
 
+// // [[Rcpp::export(".cyclopsGetLogLikelihoodGradient")]]
+// NumericVector cyclopsGetLogLikelihoodGradient(SEXP inRcppCcdInterface, int index) {
+//     using namespace bsccs;
+//
+//     XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+//
+//     auto& ccd = interface->getCcd();
+//     auto& data = interface->getModelData();
+//
+//     // const auto offset = data.getHasOffsetCovariate();
+//     const int offset = 0;
+//
+//     return ccd.getLogLikelihoodGradient(index + offset);
+// }
+
+// [[Rcpp::export(".cyclopsGetLogLikelihoodHessianDiagonal")]]
+NumericVector cyclopsGetLogLikelihoodHessianDiagonal(SEXP inRcppCcdInterface, SEXP sexpCovariates) {
+    using namespace bsccs;
+
+    XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+
+    NumericVector diagonals;
+
+    if (!Rf_isNull(sexpCovariates)) {
+
+        auto& ccd = interface->getCcd();
+        auto& data = interface->getModelData();
+
+        const std::vector<double>& bitCovariates = as<std::vector<double>>(sexpCovariates);
+        ProfileVector covariates = reinterpret_cast<const std::vector<int64_t>&>(bitCovariates);
+
+        for (ProfileVector::const_iterator it = covariates.begin();
+             it != covariates.end(); ++it) {
+
+            int index = data.getColumnIndexByName(*it);
+
+            if (index == -1) {
+                std::stringstream error;
+                error << "Variable " << *it << " not found.";
+                interface->handleError(error.str());
+            } else {
+                diagonals.push_back(-ccd.getHessianDiagonal(index));
+            }
+        }
+    }
+
+    return diagonals;
+}
+
 namespace bsccs {
 
 void RcppCcdInterface::appendRList(Rcpp::List& list, const Rcpp::List& append) {

diff --git a/src/RcppCyclopsInterface.h b/src/RcppCyclopsInterface.h
@@ -110,11 +110,10 @@ class RcppCcdInterface : public CcdInterface {
     static NoiseLevels parseNoiseLevel(const std::string& noiseName);
   	static SelectorType parseSelectorType(const std::string& selectorName);
   	static NormalizationType parseNormalizationType(const std::string& normalizationName);
+  	static void handleError(const std::string& str);
 
 protected:
 
-		static void handleError(const std::string& str);
-
     priors::JointPriorPtr makePrior(
             const std::vector<std::string>& priorName,
             bsccs::priors::PriorFunctionPtr& priorFunctionPtr,

diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -404,6 +404,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cyclopsGetLogLikelihoodHessianDiagonal
+NumericVector cyclopsGetLogLikelihoodHessianDiagonal(SEXP inRcppCcdInterface, SEXP sexpCovariates);
+RcppExport SEXP _Cyclops_cyclopsGetLogLikelihoodHessianDiagonal(SEXP inRcppCcdInterfaceSEXP, SEXP sexpCovariatesSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP);
+    Rcpp::traits::input_parameter< SEXP >::type sexpCovariates(sexpCovariatesSEXP);
+    rcpp_result_gen = Rcpp::wrap(cyclopsGetLogLikelihoodHessianDiagonal(inRcppCcdInterface, sexpCovariates));
+    return rcpp_result_gen;
+END_RCPP
+}
 // listGPUDevices
 Rcpp::CharacterVector listGPUDevices();
 RcppExport SEXP _Cyclops_listGPUDevices() {
@@ -902,6 +914,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_Cyclops_cyclopsRunBootstrap", (DL_FUNC) &_Cyclops_cyclopsRunBootstrap, 4},
     {"_Cyclops_cyclopsLogModel", (DL_FUNC) &_Cyclops_cyclopsLogModel, 1},
     {"_Cyclops_cyclopsInitializeModel", (DL_FUNC) &_Cyclops_cyclopsInitializeModel, 4},
+    {"_Cyclops_cyclopsGetLogLikelihoodHessianDiagonal", (DL_FUNC) &_Cyclops_cyclopsGetLogLikelihoodHessianDiagonal, 2},
     {"_Cyclops_listGPUDevices", (DL_FUNC) &_Cyclops_listGPUDevices, 0},
     {"_Cyclops_getDefaultGPUDevice", (DL_FUNC) &_Cyclops_getDefaultGPUDevice, 0},
     {"_Cyclops_isSorted", (DL_FUNC) &_Cyclops_isSorted, 3},

diff --git a/tests/testthat/test-finiteMLE.R b/tests/testthat/test-finiteMLE.R
@@ -25,7 +25,8 @@ test_that("Check for infinite MLE in Cox example with no outcomes in one treatme
     cyclopsData <- createCyclopsData(Surv(time, outcome) ~ exposure, data = data, modelType = "cox")
     expect_warning(fit <- fitCyclopsModel(cyclopsData), regexp =".*coefficient may be infinite.*")
 
-    ci <- confint(fit, parm = "exposureTRUE", level = 0.9)
+    expect_error(ci <- confint(fit, parm = "exposureTRUE", level = 0.9))
+    ci <- confint(fit, parm = "exposureTRUE", level = 0.9, maximumCurvature = 0)
     expect_true(is.na(ci[2]))
 })