diff --git a/NEWS.md b/NEWS.md index ded39d8f..e5eded57 100644 --- 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 index 507e596d..e244e838 100644 --- 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 index e7e46c23..2fb75310 100644 --- 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 index 092427df..53d29733 100644 --- a/man/confint.cyclopsFit.Rd +++ b/man/confint.cyclopsFit.Rd @@ -11,6 +11,7 @@ overrideNoRegularization = FALSE, includePenalty = TRUE, rescale = FALSE, + maximumCurvature = -0.001, ... ) } @@ -28,6 +29,8 @@ either a vector of numbers of covariateId names} \item{rescale}{Boolean: rescale coefficients for unnormalized covariate values} +\item{maximumCurvature}{Numeric: maximum curvature of log-likelihood profile required to compute CI} + \item{...}{Additional argument(s) for methods} } \value{ diff --git a/man/cyclops.Rd b/man/cyclops.Rd index 2936c6f2..0a3c97c0 100644 --- a/man/cyclops.Rd +++ b/man/cyclops.Rd @@ -2,7 +2,8 @@ % Please edit documentation in R/cyclops.R \docType{package} \name{cyclops} -\alias{cyclops} +\alias{Cyclops} +\alias{Cyclops-package} \title{Cyclops: Cyclic coordinate descent for logistic, Poisson and survival analysis} \description{ The Cyclops package incorporates cyclic coordinate descent and diff --git a/src/RcppCyclopsInterface.cpp b/src/RcppCyclopsInterface.cpp index b8b18072..28b4044a 100644 --- 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 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 interface(inRcppCcdInterface); + + NumericVector diagonals; + + if (!Rf_isNull(sexpCovariates)) { + + auto& ccd = interface->getCcd(); + auto& data = interface->getModelData(); + + const std::vector& bitCovariates = as>(sexpCovariates); + ProfileVector covariates = reinterpret_cast&>(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 index 83808f7b..182e9ad4 100644 --- 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& priorName, bsccs::priors::PriorFunctionPtr& priorFunctionPtr, diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 260cc507..fb02ef47 100644 --- 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 index ab76f8d5..14a8942f 100644 --- 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])) })