prelim impl Schoenfeld residuals

NAMESPACE

@@ -8,6 +8,7 @@ S3method(logLik,cyclopsFit)
 S3method(predict,cyclopsFit)
 S3method(print,cyclopsData)
 S3method(print,cyclopsFit)
+S3method(residuals,cyclopsFit)
 S3method(summary,cyclopsData)
 S3method(survfit,cyclopsFit)
 S3method(vcov,cyclopsFit)
@@ -70,6 +71,7 @@ importFrom(stats,predict)
 importFrom(stats,qchisq)
 importFrom(stats,qnorm)
 importFrom(stats,rbinom)
+importFrom(stats,residuals)
 importFrom(stats,rexp)
 importFrom(stats,rnorm)
 importFrom(stats,rpois)

NEWS.md

@@ -1,3 +1,8 @@
+develop
+==============
+
+1. add Schoenfeld residual output for Cox models
+
 Cyclops v3.5.0
 ==============
 

R/ModelFit.R

@@ -897,7 +897,7 @@ confint.cyclopsFit <- function(object, parm, level = 0.95, #control,
 #'
 #' @param object    Fitted Cyclops model object
 #' @param parm      Specification of which parameter requires profiling,
-#'                  either a vector of numbers of covariateId names
+#'                  either a vector of numbers or covariateId names
 #' @param x         Vector of values of the parameter
 #' @param bounds    Pair of values to bound adaptive profiling
 #' @param tolerance Absolute tolerance allowed for adaptive profiling

R/RcppExports.R

@@ -65,6 +65,10 @@
     invisible(.Call(`_Cyclops_cyclopsLogResult`, inRcppCcdInterface, fileName, withASE))
 }
 
+.cyclopsGetSchoenfeldResiduals <- function(inRcppCcdInterface, sexpBitCovariates) {
+    .Call(`_Cyclops_cyclopsGetSchoenfeldResiduals`, inRcppCcdInterface, sexpBitCovariates)
+}
+
 .cyclopsGetFisherInformation <- function(inRcppCcdInterface, sexpBitCovariates) {
     .Call(`_Cyclops_cyclopsGetFisherInformation`, inRcppCcdInterface, sexpBitCovariates)
 }

R/Residuals.R

@@ -0,0 +1,51 @@
+# @file Residuals.R
+#
+# Copyright 2024 Observational Health Data Sciences and Informatics
+#
+# This file is part of cyclops
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+#' @method residuals cyclopsFit
+#' @title Model residuals
+#'
+#' @description
+#' \code{residuals.cyclopsFit} computes model residuals for Cox model-based Cyclops objects
+#'
+#' @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 type      Character string indicating the type of residual desires. Possible
+#'                  values are "schoenfeld".
+#'
+#' @importFrom stats residuals
+#'
+#' @export
+residuals.cyclopsFit <- function(object, parm, type = "schoenfeld", ...) {
+    modelType <- object$cyclopsData$modelType
+    if (modelType != "cox") {
+        stop("Residuals for only Cox models are implemented")
+    }
+    if (type != "schoenfeld") {
+        stop("Only Schoenfeld residuals are implemented")
+    }
+
+    .checkInterface(object$cyclopsData, testOnly = TRUE)
+
+    res <- .cyclopsGetSchoenfeldResiduals(cyclopsFitRight$interface, NULL)
+
+    result <- res$residuals
+    names(result) <- res$times
+
+    return(rev(result))
+}

src/RcppCyclopsInterface.cpp

@@ -214,6 +214,39 @@ 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<IdType> indices;
+    if (!Rf_isNull(sexpBitCovariates)) {
+        const std::vector<double>& bitCovariates = as<std::vector<double>>(sexpBitCovariates);
+        ProfileVector covariates =  reinterpret_cast<const std::vector<int64_t>&>(bitCovariates); // as<ProfileVector>(sexpCovariates);
+        for (auto it = covariates.begin(); it != covariates.end(); ++it) {
+            size_t index = interface->getModelData().getColumnIndex(*it);
+            indices.push_back(index);
+        }
+    } else {
+        indices.push_back(0);
+    }
+    if (indices.size() != 1) {
+        Rcpp::stop("Not yet implemented");
+    }
+
+    std::vector<double> residuals;
+    std::vector<double> times;
+
+    interface->getCcd().getSchoenfeldResiduals(indices[0],
+                      residuals, times);
+
+    return DataFrame::create(
+        Named("residuals") = residuals,
+        Named("times") = times
+    );
+}
+
 // [[Rcpp::export(".cyclopsGetFisherInformation")]]
 Eigen::MatrixXd cyclopsGetFisherInformation(SEXP inRcppCcdInterface, const SEXP sexpBitCovariates) {
 	using namespace bsccs;

src/RcppExports.cpp

@@ -189,6 +189,18 @@ BEGIN_RCPP
     return R_NilValue;
 END_RCPP
 }
+// cyclopsGetSchoenfeldResiduals
+Rcpp::DataFrame cyclopsGetSchoenfeldResiduals(SEXP inRcppCcdInterface, const SEXP sexpBitCovariates);
+RcppExport SEXP _Cyclops_cyclopsGetSchoenfeldResiduals(SEXP inRcppCcdInterfaceSEXP, SEXP sexpBitCovariatesSEXP) {
+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_result_gen = Rcpp::wrap(cyclopsGetSchoenfeldResiduals(inRcppCcdInterface, sexpBitCovariates));
+    return rcpp_result_gen;
+END_RCPP
+}
 // cyclopsGetFisherInformation
 Eigen::MatrixXd cyclopsGetFisherInformation(SEXP inRcppCcdInterface, const SEXP sexpBitCovariates);
 RcppExport SEXP _Cyclops_cyclopsGetFisherInformation(SEXP inRcppCcdInterfaceSEXP, SEXP sexpBitCovariatesSEXP) {
@@ -862,6 +874,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_cyclopsGetSchoenfeldResiduals", (DL_FUNC) &_Cyclops_cyclopsGetSchoenfeldResiduals, 2},
     {"_Cyclops_cyclopsGetFisherInformation", (DL_FUNC) &_Cyclops_cyclopsGetFisherInformation, 2},
     {"_Cyclops_cyclopsSetPrior", (DL_FUNC) &_Cyclops_cyclopsSetPrior, 6},
     {"_Cyclops_cyclopsTestParameterizedPrior", (DL_FUNC) &_Cyclops_cyclopsTestParameterizedPrior, 4},

src/cyclops/CyclicCoordinateDescent.cpp

@@ -1084,8 +1084,8 @@ void CyclicCoordinateDescent::findMode(
 		modelSpecifics.setPriorParams(temp);
 		//modelSpecifics.resetBeta();
 	}
-	
-	auto cycle = [this,&lastObjFunc,&lastObjFuncVec,&iteration,algorithmType,&allDelta,&doItAll] {
+
+	auto cycle = [this,&iteration,algorithmType,&allDelta,&doItAll] {
 /*
 		if (iteration%10==0) {
 			std::cout<<"iteration " << iteration << " ";
@@ -1172,7 +1172,7 @@ void CyclicCoordinateDescent::findMode(
 						if (delta != 0.0) {
 							sufficientStatisticsKnown = false;
 							updateSufficientStatistics(delta, index);
-						}	
+						}
 				}
 	    			log(index);
 	    		}
@@ -1815,6 +1815,17 @@ void CyclicCoordinateDescent::turnOffSyncCV() {
 	modelSpecifics.turnOffSyncCV();
 }
 
+void CyclicCoordinateDescent::getSchoenfeldResiduals(const IdType index,
+                                                     std::vector<double>& residuals,
+                                                     std::vector<double>& times) {
+
+    checkAllLazyFlags();
+
+    modelSpecifics.computeSchoenfeldResiduals(index,
+                                              residuals, times,
+                                              false);
+}
+
 std::vector<double> CyclicCoordinateDescent::getPredictiveLogLikelihood(std::vector<std::vector<double>>& weightsPool) {
 	xBetaKnown = false;
 	if (usingGPU && syncCV) xBetaKnown = true;

src/cyclops/CyclicCoordinateDescent.h

@@ -130,13 +130,17 @@ class CyclicCoordinateDescent {
 
 	std::vector<double> getCensorWeights(); // ESK:
 
+	void getSchoenfeldResiduals(const IdType index,
+                             std::vector<double>& residuals,
+                             std::vector<double>& times);
+
 	void setLogisticRegression(bool idoLR);
 
 //	template <typename T>
 	void setBeta(const std::vector<double>& beta);
 
 	void setStartingBeta(const std::vector<double>& inStartingBeta);
-	
+
 	void setBeta(int i, double beta);
 
 //	void double getHessianComponent(int i, int j);

src/cyclops/engine/AbstractModelSpecifics.h

@@ -77,6 +77,14 @@ class AbstractModelSpecifics {
 	virtual void computeFisherInformation(int indexOne, int indexTwo,
 			double *oinfo, bool useWeights) = 0; // pure virtual
 
+	virtual void computeSchoenfeldResiduals(int indexOne,
+                                         std::vector<double>& residuals,
+                                         std::vector<double>& times,
+                                            // double* residuals,
+                                            // double* numerators,
+                                            // double* denominators,
+                                            bool useWeights) = 0; // pure virtual
+
 	virtual void updateXBeta(double realDelta, int index, bool useWeights) = 0; // pure virtual
 
 	virtual void computeXBeta(double* beta, bool useWeights) = 0; // pure virtual

src/cyclops/engine/ModelSpecifics.h

@@ -197,6 +197,11 @@ 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,
+                                 // double* residuals, double* numerators, double* denominators,
+                                 bool useWeights);
+
 	void updateXBeta(double delta, int index, bool useWeights);
 
 	void computeRemainingStatistics(bool useWeights);
@@ -322,6 +327,12 @@ class ModelSpecifics : public AbstractModelSpecifics, BaseModel {
 	template <class IteratorTypeOne, class IteratorTypeTwo, class Weights>
 	void computeFisherInformationImpl(int indexOne, int indexTwo, double *oinfo, Weights w);
 
+	template <class IteratorType, class Weights>
+	void getSchoenfeldResidualsImpl(int index,
+                                 std::vector<double>& residuals,
+                                 std::vector<double>& times,
+                                 Weights w);
+
 	template<class IteratorType>
 	SparseIterator<RealType> getSubjectSpecificHessianIterator(int index);
 
@@ -1260,7 +1271,7 @@ struct CoxProportionalHazards : public Storage<RealType>, OrderedData, SurvivalP
 		return static_cast<RealType>(yi);
 	}
 
-    template <class IteratorType, class Weights>		
+    template <class IteratorType, class Weights>
     inline void incrementGradientAndHessian( // TODO Make static again?
             const IteratorType& it,
             Weights false_signature,