diff --git a/NAMESPACE b/NAMESPACE index fa86367c..9aa3f4ab 100644 --- 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 index 5df701b5..e7e46c23 100644 --- 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 index 9d7f5cb4..66faec1c 100644 --- 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 new file mode 100644 index 00000000..3949e40d --- /dev/null +++ b/man/testProportionality.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Residuals.R +\name{testProportionality} +\alias{testProportionality} +\title{Test hazard ratio proportionality assumption} +\usage{ +testProportionality(object, parm, transformedTimes) +} +\arguments{ +\item{object}{A Cyclops model fit object} + +\item{parm}{A specification of which parameters require residuals, +either a vector of numbers or covariateId names} + +\item{transformedTimes}{Vector of transformed time} +} +\description{ +\code{testProportionality} tests the hazard ratio proportionality assumption +of a Cyclops model fit object +} diff --git a/src/RcppCyclopsInterface.cpp b/src/RcppCyclopsInterface.cpp index a8a93646..e5998b44 100644 --- 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 interface(inRcppCcdInterface); +std::vector getIndices(XPtr& interface, const SEXP sexpBitCovariates) { + using namespace bsccs; std::vector indices; if (!Rf_isNull(sexpBitCovariates)) { const std::vector& bitCovariates = as>(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& covariate) { + using namespace bsccs; + XPtr interface(inRcppCcdInterface); + + std::vector indices = getIndices(interface, sexpBitCovariates); + std::vector residuals; + std::vector times; + + std::vector score; + score.resize(2); + // double score = 0.0; + + interface->getCcd().getSchoenfeldResiduals(indices[0], + &residuals, ×, &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 interface(inRcppCcdInterface); + + std::vector indices = getIndices(interface, sexpBitCovariates); std::vector residuals; std::vector times; interface->getCcd().getSchoenfeldResiduals(indices[0], - residuals, times); + &residuals, ×, nullptr, nullptr); return DataFrame::create( Named("residuals") = residuals, diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index b7de34a6..260cc507 100644 --- 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& 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& >::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 index 0507700c..18b564a7 100644 --- 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& residuals, - std::vector& times) { + std::vector* residuals, + std::vector* times, + std::vector* 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 index e277f41f..a0fb8403 100644 --- a/src/cyclops/CyclicCoordinateDescent.h +++ b/src/cyclops/CyclicCoordinateDescent.h @@ -131,8 +131,10 @@ class CyclicCoordinateDescent { std::vector getCensorWeights(); // ESK: void getSchoenfeldResiduals(const IdType index, - std::vector& residuals, - std::vector& times); + std::vector* residuals, + std::vector* times, + std::vector* covariate, + double* score); void setLogisticRegression(bool idoLR); diff --git a/src/cyclops/engine/AbstractModelSpecifics.h b/src/cyclops/engine/AbstractModelSpecifics.h index 020fe425..6f2425b8 100644 --- 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& residuals, - std::vector& times, + std::vector* residuals, + std::vector* 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 index 1256a50d..4b758245 100644 --- 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& residuals, std::vector& times, + std::vector* residuals, std::vector* times, + double* covariate, double* score, // double* residuals, double* numerators, double* denominators, bool useWeights); @@ -329,8 +330,10 @@ class ModelSpecifics : public AbstractModelSpecifics, BaseModel { template void getSchoenfeldResidualsImpl(int index, - std::vector& residuals, - std::vector& times, + std::vector* residuals, + std::vector* times, + double* covariate, + double* score, Weights w); template diff --git a/src/cyclops/engine/ModelSpecifics.hpp b/src/cyclops/engine/ModelSpecifics.hpp index 71b671e5..76ff940c 100644 --- a/src/cyclops/engine/ModelSpecifics.hpp +++ b/src/cyclops/engine/ModelSpecifics.hpp @@ -734,25 +734,39 @@ void ModelSpecifics::getPredictiveEstimates(double* y, doubl template template void ModelSpecifics::getSchoenfeldResidualsImpl(int index, - std::vector& residuals, - std::vector& times, + std::vector* residuals, + std::vector* times, + double* covariate, + double* score, Weights w) { - residuals.clear(); - times.clear(); - // TODO: only written for accummulive models (Cox, Fine/Grey) + const bool hasResiduals = residuals != nullptr; + const bool hasTimes = times != nullptr; + const bool hasScore = covariate != nullptr && score != nullptr; + + if (hasResiduals) { + residuals->clear(); + } + if (hasTimes) { + times->clear(); + } - // std::cerr << "start " << index << "\n"; + RealType gradient = static_cast(0); + RealType hessian = static_cast(0); - // int outcomeIndex = 0; + // TODO: only written for accummulive models (Cox, Fine/Grey) // if (sparseIndices[index] == nullptr || sparseIndices[index]->size() > 0) { // IteratorType it(sparseIndices[index].get(), N); IteratorType it(hX, index); - RealType accNumerator = static_cast(0); - RealType accDenominator = static_cast(0); + RealType resNumerator = static_cast(0); + RealType resDenominator = static_cast(0); + + RealType scoreNumerator1 = static_cast(0); + RealType scoreNumerator2 = static_cast(0); + RealType scoreDenominator = static_cast(0); // find start relavent accumulator reset point auto reset = begin(accReset); @@ -762,28 +776,48 @@ void ModelSpecifics::getSchoenfeldResidualsImpl(int index, for (; it; ) { int i = it.index(); + const RealType x = it.value(); if (*reset <= i) { - accNumerator = static_cast(0); - accDenominator = static_cast(0); + resNumerator = static_cast(0); + resDenominator = static_cast(0); + + scoreNumerator1 = static_cast(0); + scoreNumerator2 = static_cast(0); + scoreDenominator = static_cast(0); + ++reset; } - const auto expXBeta = exp(hXBeta[i]); - - accNumerator += expXBeta * it.value(); - accDenominator += expXBeta; + const auto expXBeta = offsExpXBeta[i]; // std::exp(hXBeta[i]); - // std::cerr << hXBeta[i] << " " << expXBeta << " " << accDenominator << "\n"; + resNumerator += expXBeta * x; + resDenominator += expXBeta; - if (hY[i] == 1) { - residuals.push_back(it.value() - accNumerator / accDenominator); - times.push_back(hOffs[i]); - // residuals[outcomeIndex] = it.value() - accNumerator / accDenominator; - // ++outcomeIndex; + if (hY[i] == 1 && hasResiduals) { + residuals->push_back(it.value() - resNumerator / resDenominator); + if (hasTimes) { + times->push_back(hOffs[i]); + } } - // residuals[i] = it.value() - accNumerator / accDenominator; + if (hasScore) { + + const auto xt = x * covariate[i]; + const auto numerator1 = expXBeta * xt; + const auto numerator2 = expXBeta * xt * xt; + + scoreNumerator1 += numerator1; + scoreNumerator2 += numerator2; + + const auto t = scoreNumerator1 / resDenominator; + gradient += hNWeight[i] * t; + hessian += hNWeight[i] * (scoreNumerator2 / resDenominator - t * t); + + if (hY[i] == 1) { + gradient -= xt; + } + } ++it; @@ -793,31 +827,42 @@ void ModelSpecifics::getSchoenfeldResidualsImpl(int index, for (++i; i < next; ++i) { if (*reset <= i) { - accNumerator = static_cast(0); - accDenominator = static_cast(0); + resNumerator = static_cast(0); + resDenominator = static_cast(0); + + scoreNumerator1 = static_cast(0); + scoreNumerator2 = static_cast(0); + scoreDenominator = static_cast(0); + ++reset; } - const auto expXBeta = exp(hXBeta[i]); + const auto expXBeta = offsExpXBeta[i]; // std::exp(hXBeta[i]); - accDenominator += expXBeta; + resDenominator += expXBeta; - if (hY[i] == 1) { - residuals.push_back(-accNumerator / accDenominator); - times.push_back(hOffs[i]); - // residuals[outcomeIndex] = -accNumerator / accDenominator; - // ++outcomeIndex; + if (hY[i] == 1 && hasResiduals) { + residuals->push_back(-resNumerator / resDenominator); + if (hasTimes) { + times->push_back(hOffs[i]); + } } - // residuals[i] = -accNumerator / accDenominator; + if (hasScore) { + // TODO incr gradient / Hessian + Rcpp::stop("Not yet implemented"); + } } } } - // } else { - // // TODO: fill residuals with 0s - // } - // std::cerr << "stop\n"; + if (hasScore) { + score[0] = static_cast(gradient); + score[1] = static_cast(hessian); + // *score = static_cast(gradient / hessian); + std::cerr << "score = " << gradient << " / " << hessian << "\n"; + // std::cerr << "hess2 = " << hess2 << " or " << static_cast(2.0) * hXjX[index] << "\n"; + } } // TODO The following function is an example of a double-dispatch, rewrite without need for virtual function @@ -1301,6 +1346,8 @@ void ModelSpecifics::computeGradientAndHessianImpl(int index hessian += static_cast(2.0) * hXjX[index]; } + std::cerr << "hessian = " << hessian << " @ " << index << "\n"; + *ogradient = static_cast(gradient); *ohessian = static_cast(hessian); @@ -1456,24 +1503,26 @@ void ModelSpecifics::computeThirdDerivativeImpl(int index, d template void ModelSpecifics::computeSchoenfeldResiduals(int indexOne, - std::vector& residuals, - std::vector& times, + std::vector* residuals, + std::vector* otimes, + double* covariate, + double* score, bool useWeights) { if (useWeights) { throw new std::logic_error("Weights are not yet implemented in Schoenfeld residual calculations"); } else { // no weights switch (hX.getFormatType(indexOne)) { case INDICATOR : - getSchoenfeldResidualsImpl>(indexOne, residuals, times, weighted); + getSchoenfeldResidualsImpl>(indexOne, residuals, otimes, covariate, score, weighted); break; case SPARSE : - getSchoenfeldResidualsImpl>(indexOne, residuals, times, weighted); + getSchoenfeldResidualsImpl>(indexOne, residuals, otimes, covariate, score, weighted); break; case DENSE : - getSchoenfeldResidualsImpl>(indexOne, residuals, times, weighted); + getSchoenfeldResidualsImpl>(indexOne, residuals, otimes, covariate, score, weighted); break; case INTERCEPT : - getSchoenfeldResidualsImpl>(indexOne, residuals, times, weighted); + getSchoenfeldResidualsImpl>(indexOne, residuals, otimes, covariate, score, weighted); break; } } diff --git a/tests/testthat/test-residuals.R b/tests/testthat/test-residuals.R new file mode 100644 index 00000000..d7043645 --- /dev/null +++ b/tests/testthat/test-residuals.R @@ -0,0 +1,26 @@ +library("testthat") +library("survival") + +suppressWarnings(RNGversion("3.5.0")) + +test_that("Check Schoenfeld residual tests", { + gfit <- coxph(Surv(futime, fustat) ~ age + ecog.ps, + data=ovarian) + gres <- residuals(gfit, "schoenfeld")[,1] + + data <- createCyclopsData(Surv(futime, fustat) ~ age + ecog.ps, + data=ovarian, modelType = "cox") + cfit <- fitCyclopsModel(data) + cres <- residuals(cfit, "schoenfeld") + + expect_equal(cres, gres) + + gold <- cox.zph(gfit, transform = "identity") + + summary(lm(cres ~ as.numeric(names(cres)) - 1)) + logLik(lm(cres ~ as.numeric(names(cres)) - 1)) + + testProportionality(cfit, parm = NULL, transformedTimes = ovarian$futime) + +}) +