intermediate work on score test for proportionality

OHDSI · Nov 5, 2024 · 70c5e64 · 70c5e64
1 parent cd308e6
commit 70c5e64
Show file tree

Hide file tree

Showing 12 changed files with 253 additions and 57 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -47,6 +47,7 @@ export(runBootstrap)
 export(setOpenCLDevice)
 export(simulateCyclopsData)
 export(splitTime)
+export(testProportionality)
 import(Matrix)
 import(Rcpp)
 import(dplyr)

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -65,6 +65,10 @@
     invisible(.Call(`_Cyclops_cyclopsLogResult`, inRcppCcdInterface, fileName, withASE))
 }
 
+.cyclopsTestProportionality <- function(inRcppCcdInterface, sexpBitCovariates, covariate) {
+    .Call(`_Cyclops_cyclopsTestProportionality`, inRcppCcdInterface, sexpBitCovariates, covariate)
+}
+
 .cyclopsGetSchoenfeldResiduals <- function(inRcppCcdInterface, sexpBitCovariates) {
     .Call(`_Cyclops_cyclopsGetSchoenfeldResiduals`, inRcppCcdInterface, sexpBitCovariates)
 }

diff --git a/R/Residuals.R b/R/Residuals.R
@@ -50,3 +50,35 @@ residuals.cyclopsFit <- function(object, parm, type = "schoenfeld", ...) {
 
     return(rev(result))
 }
+
+#' @title Test hazard ratio proportionality assumption
+#'
+#' @description
+#' \code{testProportionality} tests the hazard ratio proportionality assumption
+#' of a Cyclops model fit object
+#'
+#' @param object    A Cyclops model fit object
+#' @param parm      A specification of which parameters require residuals,
+#'                  either a vector of numbers or covariateId names
+#' @param transformedTimes Vector of transformed time
+#'
+#' @export
+testProportionality <- function(object, parm, transformedTimes) {
+
+    .checkInterface(object$cyclopsData, testOnly = TRUE)
+
+    if (object$cyclopsData$modelType != "cox") {
+        stop("Proportionality test for only Cox models are implemented")
+    }
+
+    if (getNumberOfRows(object$cyclopsData) != length(transformedTimes)) {
+        stop("Incorrect 'transformedTime' length")
+    }
+
+    transformedTimes <- transformedTimes[object$cyclopsData$sortOrder]
+    message("TODO: permute transformedTimes")
+
+    res <- .cyclopsTestProportionality(object$interface, NULL, transformedTimes)
+
+    return(res)
+}
diff --git a/man/testProportionality.Rd b/man/testProportionality.Rd
diff --git a/src/RcppCyclopsInterface.cpp b/src/RcppCyclopsInterface.cpp
@@ -214,12 +214,9 @@ void cyclopsLogResult(SEXP inRcppCcdInterface, const std::string& fileName, bool
     interface->logResultsToFile(fileName, withASE);
 }
 
-// [[Rcpp::export(".cyclopsGetSchoenfeldResiduals")]]
-Rcpp::DataFrame cyclopsGetSchoenfeldResiduals(SEXP inRcppCcdInterface,
-                                                  const SEXP sexpBitCovariates) {
-    using namespace bsccs;
-    XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+std::vector<bsccs::IdType> getIndices(XPtr<bsccs::RcppCcdInterface>& interface, const SEXP sexpBitCovariates) {
 
+    using namespace bsccs;
     std::vector<IdType> indices;
     if (!Rf_isNull(sexpBitCovariates)) {
         const std::vector<double>& bitCovariates = as<std::vector<double>>(sexpBitCovariates);
@@ -235,11 +232,47 @@ Rcpp::DataFrame cyclopsGetSchoenfeldResiduals(SEXP inRcppCcdInterface,
         Rcpp::stop("Not yet implemented");
     }
 
+    return indices;
+}
+
+// [[Rcpp::export(".cyclopsTestProportionality")]]
+Rcpp::List cyclopsTestProportionality(SEXP inRcppCcdInterface,
+                                      const SEXP sexpBitCovariates,
+                                      std::vector<double>& covariate) {
+    using namespace bsccs;
+    XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+
+    std::vector<IdType> indices = getIndices(interface, sexpBitCovariates);
+    std::vector<double> residuals;
+    std::vector<double> times;
+
+    std::vector<double> score;
+    score.resize(2);
+    // double score = 0.0;
+
+    interface->getCcd().getSchoenfeldResiduals(indices[0],
+                      &residuals, &times, &covariate, &score[0]);
+
+    return List::create(
+        Named("transformed") = covariate,
+        Named("score") = score,
+        Named("residuals") = residuals,
+        Named("times") = times
+    );
+}
+
+// [[Rcpp::export(".cyclopsGetSchoenfeldResiduals")]]
+Rcpp::DataFrame cyclopsGetSchoenfeldResiduals(SEXP inRcppCcdInterface,
+                                              const SEXP sexpBitCovariates) {
+    using namespace bsccs;
+    XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+
+    std::vector<IdType> indices = getIndices(interface, sexpBitCovariates);
     std::vector<double> residuals;
     std::vector<double> times;
 
     interface->getCcd().getSchoenfeldResiduals(indices[0],
-                      residuals, times);
+                      &residuals, &times, nullptr, nullptr);
 
     return DataFrame::create(
         Named("residuals") = residuals,

diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -189,6 +189,19 @@ BEGIN_RCPP
     return R_NilValue;
 END_RCPP
 }
+// cyclopsTestProportionality
+Rcpp::List cyclopsTestProportionality(SEXP inRcppCcdInterface, const SEXP sexpBitCovariates, std::vector<double>& covariate);
+RcppExport SEXP _Cyclops_cyclopsTestProportionality(SEXP inRcppCcdInterfaceSEXP, SEXP sexpBitCovariatesSEXP, SEXP covariateSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP);
+    Rcpp::traits::input_parameter< const SEXP >::type sexpBitCovariates(sexpBitCovariatesSEXP);
+    Rcpp::traits::input_parameter< std::vector<double>& >::type covariate(covariateSEXP);
+    rcpp_result_gen = Rcpp::wrap(cyclopsTestProportionality(inRcppCcdInterface, sexpBitCovariates, covariate));
+    return rcpp_result_gen;
+END_RCPP
+}
 // cyclopsGetSchoenfeldResiduals
 Rcpp::DataFrame cyclopsGetSchoenfeldResiduals(SEXP inRcppCcdInterface, const SEXP sexpBitCovariates);
 RcppExport SEXP _Cyclops_cyclopsGetSchoenfeldResiduals(SEXP inRcppCcdInterfaceSEXP, SEXP sexpBitCovariatesSEXP) {
@@ -874,6 +887,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_Cyclops_cyclopsGetNewPredictiveLogLikelihood", (DL_FUNC) &_Cyclops_cyclopsGetNewPredictiveLogLikelihood, 2},
     {"_Cyclops_cyclopsGetLogLikelihood", (DL_FUNC) &_Cyclops_cyclopsGetLogLikelihood, 1},
     {"_Cyclops_cyclopsLogResult", (DL_FUNC) &_Cyclops_cyclopsLogResult, 3},
+    {"_Cyclops_cyclopsTestProportionality", (DL_FUNC) &_Cyclops_cyclopsTestProportionality, 3},
     {"_Cyclops_cyclopsGetSchoenfeldResiduals", (DL_FUNC) &_Cyclops_cyclopsGetSchoenfeldResiduals, 2},
     {"_Cyclops_cyclopsGetFisherInformation", (DL_FUNC) &_Cyclops_cyclopsGetFisherInformation, 2},
     {"_Cyclops_cyclopsSetPrior", (DL_FUNC) &_Cyclops_cyclopsSetPrior, 6},

diff --git a/src/cyclops/CyclicCoordinateDescent.cpp b/src/cyclops/CyclicCoordinateDescent.cpp
@@ -1816,13 +1816,23 @@ void CyclicCoordinateDescent::turnOffSyncCV() {
 }
 
 void CyclicCoordinateDescent::getSchoenfeldResiduals(const IdType index,
-                                                     std::vector<double>& residuals,
-                                                     std::vector<double>& times) {
+                                                     std::vector<double>* residuals,
+                                                     std::vector<double>* times,
+                                                     std::vector<double>* covariate,
+                                                     double* score
+                                                     ) {
 
     checkAllLazyFlags();
 
+    // double* ptrCovariate = nullptr;
+    // if (covariate != nullptr) {
+    //     ptrCovariate = &(*covariate)[0];
+    // }
+
     modelSpecifics.computeSchoenfeldResiduals(index,
                                               residuals, times,
+                                              covariate != nullptr ? &(*covariate)[0] : nullptr,
+                                              score,
                                               false);
 }
 

diff --git a/src/cyclops/CyclicCoordinateDescent.h b/src/cyclops/CyclicCoordinateDescent.h
@@ -131,8 +131,10 @@ class CyclicCoordinateDescent {
 	std::vector<double> getCensorWeights(); // ESK:
 
 	void getSchoenfeldResiduals(const IdType index,
-                             std::vector<double>& residuals,
-                             std::vector<double>& times);
+                             std::vector<double>* residuals,
+                             std::vector<double>* times,
+                             std::vector<double>* covariate,
+                             double* score);
 
 	void setLogisticRegression(bool idoLR);
 

diff --git a/src/cyclops/engine/AbstractModelSpecifics.h b/src/cyclops/engine/AbstractModelSpecifics.h
@@ -78,8 +78,10 @@ class AbstractModelSpecifics {
 			double *oinfo, bool useWeights) = 0; // pure virtual
 
 	virtual void computeSchoenfeldResiduals(int indexOne,
-                                         std::vector<double>& residuals,
-                                         std::vector<double>& times,
+                                         std::vector<double>* residuals,
+                                         std::vector<double>* times,
+                                         double* covariate,
+                                         double* score,
                                             // double* residuals,
                                             // double* numerators,
                                             // double* denominators,

diff --git a/src/cyclops/engine/ModelSpecifics.h b/src/cyclops/engine/ModelSpecifics.h
@@ -198,7 +198,8 @@ class ModelSpecifics : public AbstractModelSpecifics, BaseModel {
 	void computeFisherInformation(int indexOne, int indexTwo, double *oinfo, bool useWeights);
 
 	void computeSchoenfeldResiduals(int indexOne,
-                                 std::vector<double>& residuals, std::vector<double>& times,
+                                 std::vector<double>* residuals, std::vector<double>* times,
+                                 double* covariate, double* score,
                                  // double* residuals, double* numerators, double* denominators,
                                  bool useWeights);
 
@@ -329,8 +330,10 @@ class ModelSpecifics : public AbstractModelSpecifics, BaseModel {
 
 	template <class IteratorType, class Weights>
 	void getSchoenfeldResidualsImpl(int index,
-                                 std::vector<double>& residuals,
-                                 std::vector<double>& times,
+                                 std::vector<double>* residuals,
+                                 std::vector<double>* times,
+                                 double* covariate,
+                                 double* score,
                                  Weights w);
 
 	template<class IteratorType>