diff --git a/.Rbuildignore b/.Rbuildignore index b3ba6c21..5695b3ca 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -33,3 +33,4 @@ deprecated .github/ ^\.github$ ^CRAN-RELEASE$ +.covrignore diff --git a/.covrignore b/.covrignore new file mode 100644 index 00000000..57af5369 --- /dev/null +++ b/.covrignore @@ -0,0 +1,10 @@ +src/cyclops/io/*InputReader.* +src/cyclops/drivers/Bootstrap* +src/cyclops/drivers/Hierarchy* +src/cyclops/drivers/GridSearch* +src/cyclops/drivers/Proportion* +src/cyclops/drivers/CrossValidationSelector.h +src/cyclops/drivers/AbstractSelector.h +R/cyclops.R +R/PackageMaintenance.R + diff --git a/.gitignore b/.gitignore index 6bdb4c97..66ec7364 100644 --- a/.gitignore +++ b/.gitignore @@ -30,3 +30,5 @@ extras CRAN-RELEASE .Rdata .DS_Store +*.gcno +*.gcda diff --git a/DESCRIPTION b/DESCRIPTION index 45a4d4d3..5427b294 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: Cyclops Type: Package Title: Cyclic Coordinate Descent for Logistic, Poisson and Survival Analysis -Version: 3.1.2 +Version: 3.1.2.9999 Authors@R: c( person("Marc A.", "Suchard", email = "msuchard@ucla.edu", role = c("aut","cre")), person("Martijn J.", "Schuemie", role = "aut"), @@ -44,6 +44,7 @@ LinkingTo: Rcpp, RcppEigen (>= 0.3.2) Suggests: testthat, + readr, MASS, gnm, ggplot2, diff --git a/NEWS.md b/NEWS.md index de4884c9..f78990b3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,10 @@ +Cyclops v3.1.3 + +Changes: + +1. fixed likelihood profiling when non-convex due to numerical instability +2. fixed parsing of 64-bit covariate IDs + Cyclops v3.1.2 ============== diff --git a/R/ModelFit.R b/R/ModelFit.R index 79f00152..62d393d2 100644 --- a/R/ModelFit.R +++ b/R/ModelFit.R @@ -32,6 +32,7 @@ #' @param returnEstimates Logical, return regression coefficient estimates in Cyclops model fit object #' @param startingCoefficients Vector of starting values for optimization #' @param fixedCoefficients Vector of booleans indicating if coefficient should be fix +#' @param warnings Logical, report regularization warnings #' @param computeDevice String: Name of compute device to employ; defaults to \code{"native"} C++ on CPU #' #' @return @@ -70,6 +71,7 @@ fitCyclopsModel <- function(cyclopsData, returnEstimates = TRUE, startingCoefficients = NULL, fixedCoefficients = NULL, + warnings = TRUE, computeDevice = "native") { # Delegate to control$setHook if exists @@ -104,7 +106,7 @@ fitCyclopsModel <- function(cyclopsData, if (!is.null(prior$setHook)) { - prior$setHook(cyclopsData) # Call-back + prior$setHook(cyclopsData, warnings) # Call-back } else { prior$exclude <- .checkCovariates(cyclopsData, prior$exclude) @@ -132,7 +134,7 @@ fitCyclopsModel <- function(cyclopsData, warn <- TRUE } } - if (warn) { + if (warn && warnings) { warning("Excluding intercept from regularization") } } @@ -167,6 +169,22 @@ fitCyclopsModel <- function(cyclopsData, } } + if (any(prior$priorType == "jeffreys")) { + if (Cyclops::getNumberOfCovariates(cyclopsData) > 1) { + stop("Jeffreys prior is currently only implemented for 1 covariate") + } + + covariate <- Cyclops::getCovariateIds(cyclopsData) + if (Cyclops::getCovariateTypes(cyclopsData, covariate) != "indicator") { + count <- reduce(cyclopsData, covariate, power = 0) + sum <- reduce(cyclopsData, covariate, power = 1) + mean <- sum / count + if (!(mean == 0.0 || mean == 1.0)) { + stop("Jeffreys prior is currently only implemented for indicator covariates") + } + } + } + .cyclopsSetPrior(cyclopsData$cyclopsInterfacePtr, prior$priorType, prior$variance, prior$exclude, graph, neighborhood) } @@ -247,7 +265,7 @@ fitCyclopsModel <- function(cyclopsData, } if (!is.null(cyclopsData$censorWeights)) { - if (cyclopsData$modelType != 'fgr') { + if (cyclopsData$modelType != 'fgr' && warnings) { warning(paste0("modelType = '", cyclopsData$modelType, "' does not use censorWeights. These weights will not be passed further.")) } if (length(cyclopsData$censorWeights) != getNumberOfRows(cyclopsData)) { @@ -273,6 +291,24 @@ fitCyclopsModel <- function(cyclopsData, fit <- .cyclopsFitModel(cyclopsData$cyclopsInterfacePtr) } + if (fit$return_flag == "POOR_BLR_STEP" && control$convergenceType == "gradient") { + + if (warnings) { + warning("BLR convergence criterion failed; coefficient may be infinite") + } + + control$convergenceType <- "lange" + return(fitCyclopsModel(cyclopsData = cyclopsData, + prior = prior, + control = control, + weights = weights, + forceNewObject = forceNewObject, + returnEstimates = returnEstimates, + startingCoefficients = startingCoefficients, + fixedCoefficients = fixedCoefficients, + computeDevice = computeDevice)) + } + if (returnEstimates) { estimates <- .cyclopsLogModel(cyclopsData$cyclopsInterfacePtr) fit <- c(fit, estimates) @@ -295,12 +331,14 @@ fitCyclopsModel <- function(cyclopsData, .checkCovariates <- function(cyclopsData, covariates) { if (!is.null(covariates)) { saved <- covariates + + indices <- NULL + if (inherits(covariates, "character")) { # Try to match names indices <- match(covariates, cyclopsData$coefficientNames) covariates <- getCovariateIds(cyclopsData)[indices] } - # covariates = as.numeric(covariates) if (!bit64::is.integer64(covariates)) { covariates <- bit64::as.integer64(covariates) @@ -309,6 +347,8 @@ fitCyclopsModel <- function(cyclopsData, if (any(is.na(covariates))) { stop("Unable to match all covariates: ", paste(saved, collapse = ", ")) } + + attr(covariates, "indices") <- indices } covariates } @@ -582,7 +622,7 @@ createPrior <- function(priorType, neighborhood = NULL, useCrossValidation = FALSE, forceIntercept = FALSE) { - validNames = c("none", "laplace","normal", "barupdate", "hierarchical") + validNames = c("none", "laplace","normal", "barupdate", "hierarchical", "jeffreys") stopifnot(priorType %in% validNames) if (!is.null(exclude)) { if (!inherits(exclude, "character") && @@ -660,7 +700,7 @@ getCyclopsPredictiveLogLikelihood <- function(object, weights) { } # TODO Remove code duplication with weights section of fitCyclopsModel - .cyclopsGetPredictiveLogLikelihood(object$cyclopsData$cyclopsInterfacePtr, weights) + .cyclopsGetNewPredictiveLogLikelihood(object$cyclopsData$cyclopsInterfacePtr, weights) } #' @title Get cross-validation information from a Cyclops model fit @@ -690,7 +730,7 @@ getCrossValidationInfo <- function(object) { control$seed <- as.integer(Sys.time()) } - if (is.null(control$algoritm) || is.na(control$algorithm)) { # Provide backwards compatibility + if (is.null(control$algorithm) || is.na(control$algorithm)) { # Provide backwards compatibility control$algorithm <- "ccd" } @@ -770,6 +810,7 @@ confint.cyclopsFit <- function(object, parm, level = 0.95, #control, rescale = FALSE, ...) { .checkInterface(object$cyclopsData, testOnly = TRUE) #.setControl(object$cyclopsData$cyclopsInterfacePtr, control) + parm <- .checkCovariates(object$cyclopsData, parm) if (level < 0.01 || level > 0.99) { stop("level must be between 0 and 1") @@ -787,12 +828,15 @@ confint.cyclopsFit <- function(object, parm, level = 0.95, #control, threads, threshold, overrideNoRegularization, includePenalty) + + indices <- match(parm, getCovariateIds(object$cyclopsData)) + if (!is.null(object$scale) && rescale) { - prof$lower <- prof$lower * object$scale[as.integer(parm)] - prof$upper <- prof$upper * object$scale[as.integer(parm)] + prof$lower <- prof$lower * object$scale[indices] + prof$upper <- prof$upper * object$scale[indices] } prof <- as.matrix(as.data.frame(prof)) - rownames(prof) <- object$coefficientNames[as.integer(parm)] + rownames(prof) <- object$coefficientNames[indices] qs <- c((1 - level) / 2, 1 - (1 - level) / 2) * 100 colnames(prof)[2:3] <- paste(sprintf("%.1f", qs), "%") diff --git a/R/NewDataConversion.R b/R/NewDataConversion.R index 32b5fb3e..67a9424c 100644 --- a/R/NewDataConversion.R +++ b/R/NewDataConversion.R @@ -309,13 +309,12 @@ convertToCyclopsData.tbl_dbi <- function(outcomes, (select(outcomes, .data$y) %>% distinct() %>% count() %>% collect() > 2)) { stop("Cox model only accepts one outcome type") } - - outcomes <- outcomes %>% - arrange(.data$stratumId, desc(.data$time), .data$y, .data$rowId) if (!"time" %in% colnames(covariates)) { covariates <- covariates %>% inner_join(select(outcomes, .data$rowId, .data$time, .data$y), by = "rowId") } + outcomes <- outcomes %>% + arrange(.data$stratumId, desc(.data$time), .data$y, .data$rowId) covariates <- covariates %>% arrange(.data$covariateId, .data$stratumId, desc(.data$time), .data$y, .data$rowId) } diff --git a/R/ParameterizedPrior.R b/R/ParameterizedPrior.R index 95c53221..b0e3dee8 100644 --- a/R/ParameterizedPrior.R +++ b/R/ParameterizedPrior.R @@ -41,7 +41,7 @@ createParameterizedPrior <- function(priorType, stop("Cannot perform cross validation with a flat prior") } - setHook <- function(cyclopsData) { + setHook <- function(cyclopsData, warnings) { # closure to capture arguments if (length(priorType) > 1) { if (length(priorType) != getNumberOfCovariates(cyclopsData)) { @@ -51,7 +51,9 @@ createParameterizedPrior <- function(priorType, if (priorType[1] != "none" && .cyclopsGetHasIntercept(cyclopsData) && !forceIntercept) { priorType[1] <- "none" - warning("Excluding intercept from regularization") + if (warnings) { + warning("Excluding intercept from regularization") + } } .cyclopsSetParameterizedPrior(cyclopsData$cyclopsInterfacePtr, diff --git a/R/RcppExports.R b/R/RcppExports.R index 24c97eb0..d4d5369b 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -41,10 +41,6 @@ invisible(.Call(`_Cyclops_cyclopsSetCensorWeights`, inRcppCcdInterface, weights)) } -.cyclopsGetPredictiveLogLikelihood <- function(inRcppCcdInterface, weights) { - .Call(`_Cyclops_cyclopsGetPredictiveLogLikelihood`, inRcppCcdInterface, weights) -} - .cyclopsGetNewPredictiveLogLikelihood <- function(inRcppCcdInterface, weights) { .Call(`_Cyclops_cyclopsGetNewPredictiveLogLikelihood`, inRcppCcdInterface, weights) } @@ -53,6 +49,10 @@ .Call(`_Cyclops_cyclopsGetLogLikelihood`, inRcppCcdInterface) } +.cyclopsLogResults <- function(inRcppCcdInterface, fileName, withASE) { + invisible(.Call(`_Cyclops_cyclopsLogResult`, inRcppCcdInterface, fileName, withASE)) +} + .cyclopsGetFisherInformation <- function(inRcppCcdInterface, sexpBitCovariates) { .Call(`_Cyclops_cyclopsGetFisherInformation`, inRcppCcdInterface, sexpBitCovariates) } @@ -109,18 +109,6 @@ .Call(`_Cyclops_isSortedVectorList`, vectorList, ascending) } -#' @title Print row identifiers -#' -#' @description -#' \code{printCcdRowIds} return the row identifiers in a Cyclops data object -#' -#' @param object A Cyclops data object -#' -#' @keywords internal -printCyclopsRowIds <- function(object) { - invisible(.Call(`_Cyclops_cyclopsPrintRowIds`, object)) -} - .isRcppPtrNull <- function(x) { .Call(`_Cyclops_isRcppPtrNull`, x) } diff --git a/R/Simulation.R b/R/Simulation.R index 9129cffc..b900d351 100644 --- a/R/Simulation.R +++ b/R/Simulation.R @@ -184,7 +184,7 @@ simulateCyclopsData <- function(nstrata = 200, data$rowId <- 1:nrow(data) data <- merge(data,sim$outcomes) data <- data[order(data$stratumId,data$rowId),] - formula <- as.formula(paste(c("Surv(time,y) ~ strata(stratumId)",paste("v",1:ncovars,sep="")),collapse=" + ")) + formula <- as.formula(paste(c("Surv(time,y) ~ strata(stratumId)",paste("V",1:ncovars,sep="")),collapse=" + ")) fit <- survival::coxph(formula,data=data) if (coverage) { ci <- confint(fit) @@ -198,44 +198,44 @@ simulateCyclopsData <- function(nstrata = 200, data.frame(coef = coef(fit), lbCi95 = ci[,1], ubCi95 = ci[,2]) } -.fitUsingGnm <- function(sim,coverage=TRUE){ - if (!requireNamespace("gnm")) { - stop("gnm library required") - } - start <- Sys.time() - covariates <- sim$covariates - ncovars <- max(covariates$covariateId) - nrows <- nrow(sim$outcomes) - m <- matrix(0,nrows,ncovars) - for (i in 1:nrow(covariates)){ - m[covariates$rowId[i],covariates$covariateId[i]] <- 1 - } - data <- as.data.frame(m) - - data$rowId <- 1:nrow(data) - data <- merge(data,sim$outcomes) - data <- data[order(data$stratumId,data$rowId),] - formula <- as.formula(paste(c("y ~ v1",paste("v",2:ncovars,sep="")),collapse=" + ")) - - fit = gnm::gnm(formula, family=poisson, offset=log(time), eliminate=as.factor(data$stratumId), data = data) - #Todo: figure out how to do confidence intervals correctly - confint(fit) - fit0 = gnm::gnm(y ~ 1, family=poisson, offset=log(time), eliminate=as.factor(data$stratumId), data = data) - se <- abs(coef(fit)[[1]]/qnorm(1-pchisq(deviance(fit0)-deviance(fit),1))) - - - fit <- survival::coxph(formula,data=data) - if (coverage) { - ci <- confint(fit) - } else { - ci <- matrix(0,nrow=1,ncol=2) - } - - delta <- Sys.time() - start - writeLines(paste("Analysis took", signif(delta,3), attr(delta,"units"))) - - data.frame(coef = coef(fit), lbCi95 = ci[,1], ubCi95 = ci[,2]) -} +# .fitUsingGnm <- function(sim,coverage=TRUE){ +# if (!requireNamespace("gnm")) { +# stop("gnm library required") +# } +# start <- Sys.time() +# covariates <- sim$covariates +# ncovars <- max(covariates$covariateId) +# nrows <- nrow(sim$outcomes) +# m <- matrix(0,nrows,ncovars) +# for (i in 1:nrow(covariates)){ +# m[covariates$rowId[i],covariates$covariateId[i]] <- 1 +# } +# data <- as.data.frame(m) +# +# data$rowId <- 1:nrow(data) +# data <- merge(data,sim$outcomes) +# data <- data[order(data$stratumId,data$rowId),] +# formula <- as.formula(paste(c("y ~ v1",paste("v",2:ncovars,sep="")),collapse=" + ")) +# +# fit = gnm::gnm(formula, family=poisson, offset=log(time), eliminate=as.factor(data$stratumId), data = data) +# #Todo: figure out how to do confidence intervals correctly +# confint(fit) +# fit0 = gnm::gnm(y ~ 1, family=poisson, offset=log(time), eliminate=as.factor(data$stratumId), data = data) +# se <- abs(coef(fit)[[1]]/qnorm(1-pchisq(deviance(fit0)-deviance(fit),1))) +# +# +# fit <- survival::coxph(formula,data=data) +# if (coverage) { +# ci <- confint(fit) +# } else { +# ci <- matrix(0,nrow=1,ncol=2) +# } +# +# delta <- Sys.time() - start +# writeLines(paste("Analysis took", signif(delta,3), attr(delta,"units"))) +# +# data.frame(coef = coef(fit), lbCi95 = ci[,1], ubCi95 = ci[,2]) +# } #' @title Fit simulated data @@ -280,7 +280,7 @@ fitCyclopsSimulation <- function(sim, .fitCyclopsSimulationUsingCyclops <- function(sim, model = "logistic", - regularized = TRUE, + regularized = FALSE, coverage = TRUE, includePenalty = FALSE){ if (!regularized) @@ -299,7 +299,7 @@ fitCyclopsSimulation <- function(sim, dataPtr <- convertToCyclopsData(sim$outcomes,sim$covariates,modelType = modelType,addIntercept = !stratified) - if (regularized){ + # if (regularized){ coefCyclops <- rep(0,length(sim$effectSizes$rr)) lbCi95 <- rep(0,length(sim$effectSizes$rr)) ubCi95 <- rep(0,length(sim$effectSizes$rr)) @@ -307,7 +307,12 @@ fitCyclopsSimulation <- function(sim, if (model == "survival"){ # For some reason can't get confint twice on same dataPtr object, so recreate: dataPtr <- convertToCyclopsData(sim$outcomes,sim$covariates,modelType = modelType,addIntercept = !stratified) } - cyclopsFit <- fitCyclopsModel(dataPtr,prior = createPrior("laplace",0.1,exclude=i)) + if (regularized) { + prior <- createPrior("laplace", 0.1, exclude = i) + } else { + prior <- createPrior("none") + } + cyclopsFit <- fitCyclopsModel(dataPtr,prior = prior) coefCyclops[i] <- coef(cyclopsFit)[names(coef(cyclopsFit)) == as.character(i)] if (coverage) { if (model == "survival"){ @@ -321,23 +326,23 @@ fitCyclopsSimulation <- function(sim, } } } - } else { - cyclopsFit <- fitCyclopsModel(dataPtr, prior = createPrior("none")) - coefCyclops <- data.frame(covariateId = as.integer(names(coef(cyclopsFit))),beta=coef(cyclopsFit)) - coefCyclops <- coefCyclops[order(coefCyclops$covariateId),] - coefCyclops <- coefCyclops$beta - if (coverage) { - if (model == "survival"){ - ci <- confint(cyclopsFit,parm=i,includePenalty = includePenalty) - lbCi95 <- ci[,2] - ubCi95 <- ci[,3] - } else { - ci <- confint(cyclopsFit,parm=i,includePenalty = includePenalty) - lbCi95 <- ci[,2] - ubCi95 <- ci[,3] - } - } - } + # } else { + # cyclopsFit <- fitCyclopsModel(dataPtr, prior = createPrior("none")) + # coefCyclops <- data.frame(covariateId = as.integer(names(coef(cyclopsFit))),beta=coef(cyclopsFit)) + # coefCyclops <- coefCyclops[order(coefCyclops$covariateId),] + # coefCyclops <- coefCyclops$beta + # if (coverage) { + # if (model == "survival"){ + # ci <- confint(cyclopsFit,parm=i,includePenalty = includePenalty) + # lbCi95 <- ci[,2] + # ubCi95 <- ci[,3] + # } else { + # ci <- confint(cyclopsFit,parm=i,includePenalty = includePenalty) + # lbCi95 <- ci[,2] + # ubCi95 <- ci[,3] + # } + # } + # } delta <- Sys.time() - start writeLines(paste("Analysis took", signif(delta,3), attr(delta,"units"), "(", diff --git a/man/fitCyclopsModel.Rd b/man/fitCyclopsModel.Rd index ef5bb9dc..22381607 100644 --- a/man/fitCyclopsModel.Rd +++ b/man/fitCyclopsModel.Rd @@ -13,6 +13,7 @@ fitCyclopsModel( returnEstimates = TRUE, startingCoefficients = NULL, fixedCoefficients = NULL, + warnings = TRUE, computeDevice = "native" ) } @@ -33,6 +34,8 @@ fitCyclopsModel( \item{fixedCoefficients}{Vector of booleans indicating if coefficient should be fix} +\item{warnings}{Logical, report regularization warnings} + \item{computeDevice}{String: Name of compute device to employ; defaults to \code{"native"} C++ on CPU} } \value{ diff --git a/man/printCyclopsRowIds.Rd b/man/printCyclopsRowIds.Rd deleted file mode 100644 index f24784e9..00000000 --- a/man/printCyclopsRowIds.Rd +++ /dev/null @@ -1,15 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{printCyclopsRowIds} -\alias{printCyclopsRowIds} -\title{Print row identifiers} -\usage{ -printCyclopsRowIds(object) -} -\arguments{ -\item{object}{A Cyclops data object} -} -\description{ -\code{printCcdRowIds} return the row identifiers in a Cyclops data object -} -\keyword{internal} diff --git a/src/RcppCyclopsInterface.cpp b/src/RcppCyclopsInterface.cpp index 829d757b..ae9b1914 100644 --- a/src/RcppCyclopsInterface.cpp +++ b/src/RcppCyclopsInterface.cpp @@ -49,6 +49,10 @@ namespace bsccs { {ModelType::FINE_GRAY, "fgr"}, }; +void RcppCcdInterface::logResultsToFile(const std::string& fileName, bool withASE) { + ccd->logResults(fileName.c_str(), withASE); +} + } // namespace bsccs // [[Rcpp::export(".cyclopsGetModelTypeNames")]] @@ -145,16 +149,16 @@ void cyclopsSetCensorWeights(SEXP inRcppCcdInterface, interface->getCcd().setCensorWeights(&weights[0]); } -// [[Rcpp::export(".cyclopsGetPredictiveLogLikelihood")]] -double cyclopsGetPredictiveLogLikelihood(SEXP inRcppCcdInterface, - NumericVector& weights) { - using namespace bsccs; - XPtr interface(inRcppCcdInterface); - - // return interface->getCcd().getPredictiveLogLikelihood(&weights[0]); - Rcpp::stop("No longer implemented"); - return 0.0; -} +// // [[Rcpp::export(".cyclopsGetPredictiveLogLikelihood")]] +// double cyclopsGetPredictiveLogLikelihood(SEXP inRcppCcdInterface, +// NumericVector& weights) { +// using namespace bsccs; +// XPtr interface(inRcppCcdInterface); +// +// // return interface->getCcd().getPredictiveLogLikelihood(&weights[0]); +// Rcpp::stop("No longer implemented"); +// return 0.0; +// } // [[Rcpp::export(".cyclopsGetNewPredictiveLogLikelihood")]] double cyclopsGetNewPredictiveLogLikelihood(SEXP inRcppCcdInterface, @@ -173,6 +177,15 @@ double cyclopsGetLogLikelihood(SEXP inRcppCcdInterface) { return interface->getCcd().getLogLikelihood(); } + +// [[Rcpp::export(".cyclopsLogResults")]] +void cyclopsLogResult(SEXP inRcppCcdInterface, const std::string& fileName, bool withASE) { + using namespace bsccs; + XPtr interface(inRcppCcdInterface); + + interface->logResultsToFile(fileName, withASE); +} + // [[Rcpp::export(".cyclopsGetFisherInformation")]] Eigen::MatrixXd cyclopsGetFisherInformation(SEXP inRcppCcdInterface, const SEXP sexpBitCovariates) { using namespace bsccs; @@ -657,6 +670,8 @@ bsccs::priors::PriorType RcppCcdInterface::parsePriorType(const std::string& pri priorType = NORMAL; } else if (priorName == "barupdate") { priorType = BAR_UPDATE; + } else if (priorName == "jeffreys") { + priorType = JEFFREYS; } else { handleError("Invalid prior type."); } diff --git a/src/RcppCyclopsInterface.h b/src/RcppCyclopsInterface.h index 887e81e3..8225a2d4 100644 --- a/src/RcppCyclopsInterface.h +++ b/src/RcppCyclopsInterface.h @@ -67,6 +67,8 @@ class RcppCcdInterface : public CcdInterface { return CcdInterface::runBoostrap(ccd, modelData, savedBeta); } + void logResultsToFile(const std::string& fileName, bool withASE); + void setZeroBetaAsFixed() { CcdInterface::setZeroBetaAsFixed(ccd); } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 59d381ae..67ad4bce 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -119,18 +119,6 @@ BEGIN_RCPP return R_NilValue; END_RCPP } -// cyclopsGetPredictiveLogLikelihood -double cyclopsGetPredictiveLogLikelihood(SEXP inRcppCcdInterface, NumericVector& weights); -RcppExport SEXP _Cyclops_cyclopsGetPredictiveLogLikelihood(SEXP inRcppCcdInterfaceSEXP, SEXP weightsSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP); - Rcpp::traits::input_parameter< NumericVector& >::type weights(weightsSEXP); - rcpp_result_gen = Rcpp::wrap(cyclopsGetPredictiveLogLikelihood(inRcppCcdInterface, weights)); - return rcpp_result_gen; -END_RCPP -} // cyclopsGetNewPredictiveLogLikelihood double cyclopsGetNewPredictiveLogLikelihood(SEXP inRcppCcdInterface, NumericVector& weights); RcppExport SEXP _Cyclops_cyclopsGetNewPredictiveLogLikelihood(SEXP inRcppCcdInterfaceSEXP, SEXP weightsSEXP) { @@ -154,6 +142,18 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// cyclopsLogResult +void cyclopsLogResult(SEXP inRcppCcdInterface, const std::string& fileName, bool withASE); +RcppExport SEXP _Cyclops_cyclopsLogResult(SEXP inRcppCcdInterfaceSEXP, SEXP fileNameSEXP, SEXP withASESEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP); + Rcpp::traits::input_parameter< const std::string& >::type fileName(fileNameSEXP); + Rcpp::traits::input_parameter< bool >::type withASE(withASESEXP); + cyclopsLogResult(inRcppCcdInterface, fileName, withASE); + return R_NilValue; +END_RCPP +} // cyclopsGetFisherInformation Eigen::MatrixXd cyclopsGetFisherInformation(SEXP inRcppCcdInterface, const SEXP sexpBitCovariates); RcppExport SEXP _Cyclops_cyclopsGetFisherInformation(SEXP inRcppCcdInterfaceSEXP, SEXP sexpBitCovariatesSEXP) { @@ -353,16 +353,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// cyclopsPrintRowIds -void cyclopsPrintRowIds(Environment object); -RcppExport SEXP _Cyclops_cyclopsPrintRowIds(SEXP objectSEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Environment >::type object(objectSEXP); - cyclopsPrintRowIds(object); - return R_NilValue; -END_RCPP -} // isRcppPtrNull bool isRcppPtrNull(SEXP x); RcppExport SEXP _Cyclops_isRcppPtrNull(SEXP xSEXP) { @@ -780,9 +770,9 @@ static const R_CallMethodDef CallEntries[] = { {"_Cyclops_cyclopsGetIsRegularized", (DL_FUNC) &_Cyclops_cyclopsGetIsRegularized, 2}, {"_Cyclops_cyclopsSetWeights", (DL_FUNC) &_Cyclops_cyclopsSetWeights, 2}, {"_Cyclops_cyclopsSetCensorWeights", (DL_FUNC) &_Cyclops_cyclopsSetCensorWeights, 2}, - {"_Cyclops_cyclopsGetPredictiveLogLikelihood", (DL_FUNC) &_Cyclops_cyclopsGetPredictiveLogLikelihood, 2}, {"_Cyclops_cyclopsGetNewPredictiveLogLikelihood", (DL_FUNC) &_Cyclops_cyclopsGetNewPredictiveLogLikelihood, 2}, {"_Cyclops_cyclopsGetLogLikelihood", (DL_FUNC) &_Cyclops_cyclopsGetLogLikelihood, 1}, + {"_Cyclops_cyclopsLogResult", (DL_FUNC) &_Cyclops_cyclopsLogResult, 3}, {"_Cyclops_cyclopsGetFisherInformation", (DL_FUNC) &_Cyclops_cyclopsGetFisherInformation, 2}, {"_Cyclops_cyclopsSetPrior", (DL_FUNC) &_Cyclops_cyclopsSetPrior, 6}, {"_Cyclops_cyclopsTestParameterizedPrior", (DL_FUNC) &_Cyclops_cyclopsTestParameterizedPrior, 4}, @@ -797,7 +787,6 @@ static const R_CallMethodDef CallEntries[] = { {"_Cyclops_cyclopsInitializeModel", (DL_FUNC) &_Cyclops_cyclopsInitializeModel, 4}, {"_Cyclops_isSorted", (DL_FUNC) &_Cyclops_isSorted, 3}, {"_Cyclops_isSortedVectorList", (DL_FUNC) &_Cyclops_isSortedVectorList, 2}, - {"_Cyclops_cyclopsPrintRowIds", (DL_FUNC) &_Cyclops_cyclopsPrintRowIds, 1}, {"_Cyclops_isRcppPtrNull", (DL_FUNC) &_Cyclops_isRcppPtrNull, 1}, {"_Cyclops_cyclopsGetNumberOfStrata", (DL_FUNC) &_Cyclops_cyclopsGetNumberOfStrata, 1}, {"_Cyclops_cyclopsGetCovariateIds", (DL_FUNC) &_Cyclops_cyclopsGetCovariateIds, 1}, diff --git a/src/RcppModelData.cpp b/src/RcppModelData.cpp index c61b5898..7eb7796b 100644 --- a/src/RcppModelData.cpp +++ b/src/RcppModelData.cpp @@ -78,20 +78,20 @@ XPtr parseEnvironmentForPtr(const Environment& x) { // return ptr; // } -//' @title Print row identifiers -//' -//' @description -//' \code{printCcdRowIds} return the row identifiers in a Cyclops data object -//' -//' @param object A Cyclops data object -//' -//' @keywords internal -// [[Rcpp::export("printCyclopsRowIds")]] -void cyclopsPrintRowIds(Environment object) { - XPtr data = parseEnvironmentForPtr(object); -// std::ostreamstring stream; -// std::vector& rowsIds = data->get -} +// // ' @title Print row identifiers +// // ' +// // ' @description +// // ' \code{printCcdRowIds} return the row identifiers in a Cyclops data object +// // ' +// // ' @param object A Cyclops data object +// // ' +// // ' @keywords internal +// // [[Rcpp::export("printCyclopsRowIds")]] +// // void cyclopsPrintRowIds(Environment object) { +// // XPtr data = parseEnvironmentForPtr(object); +// // std::ostreamstring stream; +// // std::vector& rowsIds = data->get +// // } // void testCcdCode(int position) { // std::vector v{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; diff --git a/src/cyclops/CcdInterface.cpp b/src/cyclops/CcdInterface.cpp index faddefc2..5431047a 100644 --- a/src/cyclops/CcdInterface.cpp +++ b/src/cyclops/CcdInterface.cpp @@ -546,7 +546,7 @@ double CcdInterface::evaluateProfileModel(CyclicCoordinateDescent *ccd, Abstract int nThreads = (inThreads == -1) ? bsccs::thread::hardware_concurrency() : inThreads; - double mode = ccd->getLogLikelihood(); // TODO Remove + ccd->getLogLikelihood(); // TODO Remove std::ostringstream stream2; stream2 << "Using " << nThreads << " thread(s)"; diff --git a/src/cyclops/CyclicCoordinateDescent.cpp b/src/cyclops/CyclicCoordinateDescent.cpp index 7ec11e3e..fbb03df2 100644 --- a/src/cyclops/CyclicCoordinateDescent.cpp +++ b/src/cyclops/CyclicCoordinateDescent.cpp @@ -535,7 +535,7 @@ void CyclicCoordinateDescent::setCensorWeights(double* dWeights) { } double CyclicCoordinateDescent::getLogPrior(void) { - return jointPrior->logDensity(hBeta); + return jointPrior->logDensity(hBeta, *this); } double CyclicCoordinateDescent::getObjectiveFunction(int convergenceType) { @@ -844,6 +844,7 @@ bool CyclicCoordinateDescent::performCheckConvergence(int convergenceType, bool done = false; double conv; bool illconditioned = false; + bool poorBlr = false; if (convergenceType < ZHANG_OLES) { double thisObjFunc = getObjectiveFunction(convergenceType); if (thisObjFunc != thisObjFunc) { @@ -857,6 +858,16 @@ bool CyclicCoordinateDescent::performCheckConvergence(int convergenceType, illconditioned = true; } else { conv = computeConvergenceCriterion(thisObjFunc, *lastObjFunc); + if (conv == 0.0 && iteration == 1 && convergenceType != ONE_STEP && + !std::all_of(std::begin(fixBeta), std::end(fixBeta), [](bool fixed) { return fixed; }) + ) { + std::ostringstream stream; + stream << "\nWarning: BLR gradient is ill-conditioned\n" + << "Enforcing convergence!"; + logger->writeLine(stream); + illconditioned = true; + poorBlr = true; + } } *lastObjFunc = thisObjFunc; } else { // ZHANG_OLES @@ -880,7 +891,9 @@ bool CyclicCoordinateDescent::performCheckConvergence(int convergenceType, } if (epsilon > 0 && conv < epsilon) { - if (illconditioned) { + if (poorBlr) { + lastReturnFlag = POOR_BLR_STEP; + } else if (illconditioned) { lastReturnFlag = ILLCONDITIONED; } else { if (noiseLevel > SILENT) { @@ -1057,148 +1070,148 @@ void CyclicCoordinateDescent::findMode( qnQ = 0; if (qnQ > 0) { // Use quasi-Newton - using namespace Eigen; - - Eigen::MatrixXd secantsU(J, qnQ); - Eigen::MatrixXd secantsV(J, qnQ); - - VectorXd x(J); - - // Fill initial secants - int countU = 0; - int countV = 0; - - for (int q = 0; q < qnQ; ++q) { - x = Map(hBeta.data(), J); // Make copy - - cycle(); - done = check(); - if (done) break; - - if (countU == 0) { // First time through - secantsU.col(countU) = Map(hBeta.data(), J) - x; - ++countU; - } else if (countU < qnQ - 1) { // Middle case - secantsU.col(countU) = Map(hBeta.data(), J) - x; - secantsV.col(countV) = secantsU.col(countU); - ++countU; - ++countV; - } else { // Last time through - secantsV.col(countV) = Map(hBeta.data(), J) - x; - ++countV; - } - } - - // std::cerr << secantsU << std::endl << std::endl; - // std::cerr << secantsV << std::endl; - - int newestSecant = qnQ - 1; - int previousSecant = newestSecant - 1; - - while (!done) { - - // 2 cycles for each QN step - x = Map(hBeta.data(), J); // Make copy - cycle(); - // done = check(); - // if (done) { - // std::cerr << "break A" << std::endl; - // break; - // } - - secantsU.col(newestSecant) = Map(hBeta.data(), J) - x; - - VectorXd Fx = Map(hBeta.data(), J); // TODO Can remove? - - x = Map(hBeta.data(), J); - cycle(); - // done = check(); - // if (done) { - // std::cerr << "break B" << std::endl; - // break; - // } - - secantsV.col(newestSecant) = Map(hBeta.data(), J) - x; - - // Do QN step here - - // MatrixXd UtU = secantsU.transpose() * secantsU; - // MatrixXd UtV = secantsU.transpose() * secantsV; - // MatrixXd M = UtU - UtV; - - auto M = secantsU.transpose() * (secantsU - secantsV); - - auto Minv = M.inverse(); - - auto A = secantsU.transpose() * secantsU.col(newestSecant); - auto B = Minv * A; - auto C = secantsV * B; - - // MatrixXd A = secantsV * Minv; - // MatrixXd B = A * secantsU.transpose(); - // MatrixXd C = B * secantsV.col(newestSecant); // - (x - Fx) - // MatrixXd C = B * secantsU.col(newestSecant); // TODO Can remove? - - // VectorXd xqn = Map(hBeta.data(), J) + C; - VectorXd xqn = Fx + C; // TODO Can remove? - - // Save CCD solution - x = Map(hBeta.data(), J); - double ccdObjective = getLogLikelihood() + getLogPrior(); - // double savedLastObjFunc = lastObjFunc; - - Map(hBeta.data(), J) = xqn; // Set QN solution - // for (int j = 0; j < J; ++j) { - // if (sign(hBeta[j]) == xqn(j)) { - // hBeta[j] = xqn(j); - // } else { - // hBeta[j] = 0.0; - // } - // } - - // xBetaKnown = false; - // checkAllLazyFlags(); - modelSpecifics.computeXBeta(hBeta.data(), useCrossValidation); - - - // done = check(); - // if (done) { - // std::cerr << "break B" << std::endl; - // break; - // } - - double qnObjective = getLogLikelihood() + getLogPrior(); - - if (ccdObjective > qnObjective) { // Revert - Map(hBeta.data(), J) = x; // Set CCD solution - modelSpecifics.computeXBeta(hBeta.data(), useCrossValidation); - // xBetaKnown = false; - // checkAllLazyFlags(); - - double ccd2Objective = getLogLikelihood() + getLogPrior(); - - if (ccdObjective != ccd2Objective) { - std::cerr << "Poor revert: " << ccdObjective << " != " << ccd2Objective << " diff: "<< (ccdObjective - ccd2Objective) << std::endl; - } - // lastObjFunc = savedLastObjFunc; - std::cerr << "revert" << std::endl; - } else { - std::cerr << "accept" << std::endl; - // done = check(); - // if (done) { - // std::cerr << "break C" << std::endl; - // break; - // } - } - - done = check(); - - // lastObjFunc = getLogLikelihood() + getLogPrior(); - - // Get ready for next secant-pair - previousSecant = newestSecant; - newestSecant = (newestSecant + 1) % qnQ; - } + // using namespace Eigen; + // + // Eigen::MatrixXd secantsU(J, qnQ); + // Eigen::MatrixXd secantsV(J, qnQ); + // + // VectorXd x(J); + // + // // Fill initial secants + // int countU = 0; + // int countV = 0; + // + // for (int q = 0; q < qnQ; ++q) { + // x = Map(hBeta.data(), J); // Make copy + // + // cycle(); + // done = check(); + // if (done) break; + // + // if (countU == 0) { // First time through + // secantsU.col(countU) = Map(hBeta.data(), J) - x; + // ++countU; + // } else if (countU < qnQ - 1) { // Middle case + // secantsU.col(countU) = Map(hBeta.data(), J) - x; + // secantsV.col(countV) = secantsU.col(countU); + // ++countU; + // ++countV; + // } else { // Last time through + // secantsV.col(countV) = Map(hBeta.data(), J) - x; + // ++countV; + // } + // } + // + // // std::cerr << secantsU << std::endl << std::endl; + // // std::cerr << secantsV << std::endl; + // + // int newestSecant = qnQ - 1; + // int previousSecant = newestSecant - 1; + // + // while (!done) { + // + // // 2 cycles for each QN step + // x = Map(hBeta.data(), J); // Make copy + // cycle(); + // // done = check(); + // // if (done) { + // // std::cerr << "break A" << std::endl; + // // break; + // // } + // + // secantsU.col(newestSecant) = Map(hBeta.data(), J) - x; + // + // VectorXd Fx = Map(hBeta.data(), J); // TODO Can remove? + // + // x = Map(hBeta.data(), J); + // cycle(); + // // done = check(); + // // if (done) { + // // std::cerr << "break B" << std::endl; + // // break; + // // } + // + // secantsV.col(newestSecant) = Map(hBeta.data(), J) - x; + // + // // Do QN step here + // + // // MatrixXd UtU = secantsU.transpose() * secantsU; + // // MatrixXd UtV = secantsU.transpose() * secantsV; + // // MatrixXd M = UtU - UtV; + // + // auto M = secantsU.transpose() * (secantsU - secantsV); + // + // auto Minv = M.inverse(); + // + // auto A = secantsU.transpose() * secantsU.col(newestSecant); + // auto B = Minv * A; + // auto C = secantsV * B; + // + // // MatrixXd A = secantsV * Minv; + // // MatrixXd B = A * secantsU.transpose(); + // // MatrixXd C = B * secantsV.col(newestSecant); // - (x - Fx) + // // MatrixXd C = B * secantsU.col(newestSecant); // TODO Can remove? + // + // // VectorXd xqn = Map(hBeta.data(), J) + C; + // VectorXd xqn = Fx + C; // TODO Can remove? + // + // // Save CCD solution + // x = Map(hBeta.data(), J); + // double ccdObjective = getLogLikelihood() + getLogPrior(); + // // double savedLastObjFunc = lastObjFunc; + // + // Map(hBeta.data(), J) = xqn; // Set QN solution + // // for (int j = 0; j < J; ++j) { + // // if (sign(hBeta[j]) == xqn(j)) { + // // hBeta[j] = xqn(j); + // // } else { + // // hBeta[j] = 0.0; + // // } + // // } + // + // // xBetaKnown = false; + // // checkAllLazyFlags(); + // modelSpecifics.computeXBeta(hBeta.data(), useCrossValidation); + // + // + // // done = check(); + // // if (done) { + // // std::cerr << "break B" << std::endl; + // // break; + // // } + // + // double qnObjective = getLogLikelihood() + getLogPrior(); + // + // if (ccdObjective > qnObjective) { // Revert + // Map(hBeta.data(), J) = x; // Set CCD solution + // modelSpecifics.computeXBeta(hBeta.data(), useCrossValidation); + // // xBetaKnown = false; + // // checkAllLazyFlags(); + // + // double ccd2Objective = getLogLikelihood() + getLogPrior(); + // + // if (ccdObjective != ccd2Objective) { + // std::cerr << "Poor revert: " << ccdObjective << " != " << ccd2Objective << " diff: "<< (ccdObjective - ccd2Objective) << std::endl; + // } + // // lastObjFunc = savedLastObjFunc; + // std::cerr << "revert" << std::endl; + // } else { + // std::cerr << "accept" << std::endl; + // // done = check(); + // // if (done) { + // // std::cerr << "break C" << std::endl; + // // break; + // // } + // } + // + // done = check(); + // + // // lastObjFunc = getLogLikelihood() + getLogPrior(); + // + // // Get ready for next secant-pair + // previousSecant = newestSecant; + // newestSecant = (newestSecant + 1) % qnQ; + // } } else { // No QN while (!done) { cycle(); @@ -1246,6 +1259,18 @@ double CyclicCoordinateDescent::getHessianDiagonal(int index) { return g_d2; } +double CyclicCoordinateDescent::getJerkDiagonal(int index) { + + checkAllLazyFlags(); + + double jerk; + + computeNumeratorForGradient(index); // TODO Is this necessary? + modelSpecifics.computeThirdDerivative(index, &jerk, useCrossValidation); + + return jerk; +} + double CyclicCoordinateDescent::getAsymptoticVariance(int indexOne, int indexTwo) { checkAllLazyFlags(); if (!fisherInformationKnown) { @@ -1268,22 +1293,22 @@ double CyclicCoordinateDescent::getAsymptoticVariance(int indexOne, int indexTwo } } -double CyclicCoordinateDescent::getAsymptoticPrecision(int indexOne, int indexTwo) { - checkAllLazyFlags(); - if (!fisherInformationKnown) { - computeAsymptoticPrecisionMatrix(); - fisherInformationKnown = true; - } - - IndexMap::iterator itOne = hessianIndexMap.find(indexOne); - IndexMap::iterator itTwo = hessianIndexMap.find(indexTwo); - - if (itOne == hessianIndexMap.end() || itTwo == hessianIndexMap.end()) { - return NAN; - } else { - return hessianMatrix(itOne->second, itTwo->second); - } -} +// double CyclicCoordinateDescent::getAsymptoticPrecision(int indexOne, int indexTwo) { +// checkAllLazyFlags(); +// if (!fisherInformationKnown) { +// computeAsymptoticPrecisionMatrix(); +// fisherInformationKnown = true; +// } +// +// IndexMap::iterator itOne = hessianIndexMap.find(indexOne); +// IndexMap::iterator itTwo = hessianIndexMap.find(indexTwo); +// +// if (itOne == hessianIndexMap.end() || itTwo == hessianIndexMap.end()) { +// return NAN; +// } else { +// return hessianMatrix(itOne->second, itTwo->second); +// } +// } CyclicCoordinateDescent::Matrix CyclicCoordinateDescent::computeFisherInformation(const std::vector& indices) const { Matrix fisherInformation(indices.size(), indices.size()); @@ -1382,7 +1407,7 @@ void CyclicCoordinateDescent::mmUpdateAllBeta(std::vector& delta, gh[j].second /= scale; - delta[j] = jointPrior->getDelta(gh[j], hBeta, j); + delta[j] = jointPrior->getDelta(gh[j], hBeta, j, *this); // TODO this is only correct when joints are independent across dimensions } else { delta[j] = 0.0; @@ -1409,7 +1434,7 @@ double CyclicCoordinateDescent::ccdUpdateBeta(int index) { gh.second = 0.0; } - return jointPrior->getDelta(gh, hBeta, index); + return jointPrior->getDelta(gh, hBeta, index, *this); } void CyclicCoordinateDescent::axpyXBeta(const double beta, const int j) { diff --git a/src/cyclops/CyclicCoordinateDescent.h b/src/cyclops/CyclicCoordinateDescent.h index 2e4e3575..515adf21 100644 --- a/src/cyclops/CyclicCoordinateDescent.h +++ b/src/cyclops/CyclicCoordinateDescent.h @@ -97,9 +97,11 @@ class CyclicCoordinateDescent { double getHessianDiagonal(int index); + double getJerkDiagonal(int index); + double getAsymptoticVariance(int i, int j); - double getAsymptoticPrecision(int i, int j); + // double getAsymptoticPrecision(int i, int j); // void setZeroBetaFixed(void); diff --git a/src/cyclops/Types.h b/src/cyclops/Types.h index 94654bb2..32d514ac 100644 --- a/src/cyclops/Types.h +++ b/src/cyclops/Types.h @@ -102,7 +102,8 @@ enum PriorType { NONE = 0, LAPLACE, NORMAL, - BAR_UPDATE + BAR_UPDATE, + JEFFREYS }; } // namespace priors @@ -126,7 +127,8 @@ enum UpdateReturnFlags { FAIL, MAX_ITERATIONS, ILLCONDITIONED, - MISSING_COVARIATES + MISSING_COVARIATES, + POOR_BLR_STEP }; typedef std::vector ProfileVector; diff --git a/src/cyclops/engine/AbstractModelSpecifics.h b/src/cyclops/engine/AbstractModelSpecifics.h index c2c9fd7f..a0a4fb6d 100644 --- a/src/cyclops/engine/AbstractModelSpecifics.h +++ b/src/cyclops/engine/AbstractModelSpecifics.h @@ -67,6 +67,8 @@ class AbstractModelSpecifics { virtual void computeGradientAndHessian(int index, double *ogradient, double *ohessian, bool useWeights) = 0; // pure virtual + virtual void computeThirdDerivative(int index, double *othird, bool useWeights) = 0; // pure virtual + virtual void computeMMGradientAndHessian(std::vector& gh, const std::vector& fixBeta, bool useWeights) = 0; // pure virtual virtual void computeNumeratorForGradient(int index, bool useWeights) = 0; // pure virtual diff --git a/src/cyclops/engine/ModelSpecifics.h b/src/cyclops/engine/ModelSpecifics.h index 3b726a50..88dfaa68 100644 --- a/src/cyclops/engine/ModelSpecifics.h +++ b/src/cyclops/engine/ModelSpecifics.h @@ -76,6 +76,8 @@ class ModelSpecifics : public AbstractModelSpecifics, BaseModel { void computeGradientAndHessian(int index, double *ogradient, double *ohessian, bool useWeights); + void computeThirdDerivative(int index, double *othird, bool useWeights); + virtual void computeMMGradientAndHessian( std::vector& gh, const std::vector& fixBeta, @@ -244,6 +246,11 @@ class ModelSpecifics : public AbstractModelSpecifics, BaseModel { double *gradient, double *hessian, Weights w); + template + void computeThirdDerivativeImpl( + int index, + double *third, Weights w); + template void computeMMGradientAndHessianImpl( int index, double *ogradient, @@ -515,6 +522,18 @@ struct Survival { // throw new std::logic_error("Not model-specific"); } + + template + inline void incrementThirdDerivative( + const IteratorType& it, + Weights false_signature, + RealType* third, + RealType numer, RealType numer2, RealType denom, + RealType nEvents, + RealType x, RealType xBeta, RealType y) { + + throw new std::logic_error("3rd derivatives are not yet implemented"); + } }; template @@ -583,6 +602,18 @@ struct Logistic { } } + template + inline void incrementThirdDerivative( + const IteratorType& it, + Weights false_signature, + RealType* third, + RealType numer, RealType numer2, RealType denom, + RealType nEvents, + RealType x, RealType xBeta, RealType y) { + + throw new std::logic_error("3rd derivatives are not yet implemented"); + } + template static inline Fraction incrementGradientAndHessian(const Fraction& lhs, RealType numerator, RealType numerator2, RealType denominator, RealType weight, @@ -751,6 +782,18 @@ struct ConditionalPoissonRegression : public Storage, GroupedData, GLM } } + template + inline void incrementThirdDerivative( + const IteratorType& it, + Weights false_signature, + RealType* third, + RealType numer, RealType numer2, RealType denom, + RealType nEvents, + RealType x, RealType xBeta, RealType y) { + + throw new std::logic_error("3rd derivatives are not yet implemented"); + } + template static inline Fraction incrementGradientAndHessian(const Fraction& lhs, RealType numerator, RealType numerator2, RealType denominator, RealType weight, @@ -1079,7 +1122,7 @@ struct CoxProportionalHazards : public Storage, OrderedData, SurvivalP } template - void incrementGradientAndHessian( // TODO Make static again? + inline void incrementGradientAndHessian( // TODO Make static again? const IteratorType& it, Weights false_signature, RealType* gradient, RealType* hessian, @@ -1201,6 +1244,24 @@ struct BreslowTiedCoxProportionalHazards : public Storage, OrderedWith } } + template + inline void incrementThirdDerivative( + const IteratorType& it, + Weights false_signature, + RealType* third, + RealType numer, RealType numer2, RealType denom, + RealType nEvents, + RealType x, RealType xBeta, RealType y) { + + const RealType t = numer / denom; + const RealType g = nEvents * t; // Always use weights (not censured indicator) + // if (IteratorType::isIndicator) { + *third += g * (static_cast(1) - static_cast(2) * t) * (static_cast(1.0) - t); +// } else { +// throw new std::logic_error("Dense-covariate 3rd derivatives are not yet implemented"); +// } + } + template static inline Fraction incrementGradientAndHessian(const Fraction& lhs, RealType numerator, RealType numerator2, RealType denominator, RealType weight, @@ -1303,8 +1364,6 @@ struct BreslowTiedFineGray: public Storage, OrderedWithTiesData, Survi RealType numerator, RealType numerator2, RealType denominator, RealType weight, RealType xBeta, RealType y) { -// std::cout << "TODO" << std::endl; -// std::exit(-1); // breslow cox throw new std::logic_error("breslow cox model not yet support"); @@ -1427,6 +1486,17 @@ struct LeastSquares : public Storage, IndependentData, FixedPid, NoFix } } + template + inline void incrementThirdDerivative( + const IteratorType& it, + Weights false_signature, + RealType* third, + RealType numer, RealType numer2, RealType denom, + RealType nEvents, + RealType x, RealType xBeta, RealType y) { + // Do nothing + } + template inline void incrementMMGradientAndHessian( RealType& gradient, RealType& hessian, @@ -1542,6 +1612,18 @@ struct PoissonRegression : public Storage, IndependentData, GLMProject // } } + template + inline void incrementThirdDerivative( + const IteratorType& it, + Weights false_signature, + RealType* third, + RealType numer, RealType numer2, RealType denom, + RealType nEvents, + RealType x, RealType xBeta, RealType y) { + + throw new std::logic_error("3rd derivatives are not yet implemented"); + } + template inline void incrementMMGradientAndHessian( RealType& gradient, RealType& hessian, @@ -1673,9 +1755,6 @@ struct TestGradientKernel { const auto denominator = boost::get<0>(tuple); const auto weight = boost::get<1>(tuple); - // std::cerr << "N n1: " << numerator.first << " n2: " << numerator.second - // << " d: " << denominator << " w: " << weight << std::endl; - return BaseModel::template incrementGradientAndHessian( lhs, diff --git a/src/cyclops/engine/ModelSpecifics.hpp b/src/cyclops/engine/ModelSpecifics.hpp index 679e36e1..84f9eed3 100644 --- a/src/cyclops/engine/ModelSpecifics.hpp +++ b/src/cyclops/engine/ModelSpecifics.hpp @@ -1198,7 +1198,147 @@ void ModelSpecifics::computeGradientAndHessianImpl(int index #endif #endif - } +} + +template +void ModelSpecifics::computeThirdDerivative(int index, double *othird, bool useWeights) { + +#ifdef CYCLOPS_DEBUG_TIMING +#ifndef CYCLOPS_DEBUG_TIMING_LOW + auto start = bsccs::chrono::steady_clock::now(); +#endif +#endif + + if (hX.getNumberOfNonZeroEntries(index) == 0) { + *othird = 0.0; + return; + } + + // Run-time dispatch, so virtual call should not effect speed + if (useWeights) { + switch (hX.getFormatType(index)) { + case INDICATOR : + computeThirdDerivativeImpl>(index, othird, weighted); + break; + case SPARSE : + computeThirdDerivativeImpl>(index, othird, weighted); + break; + case DENSE : + computeThirdDerivativeImpl>(index, othird, weighted); + break; + case INTERCEPT : + computeThirdDerivativeImpl>(index, othird, weighted); + break; + } + } else { + switch (hX.getFormatType(index)) { + case INDICATOR : + computeThirdDerivativeImpl>(index, othird, unweighted); + break; + case SPARSE : + computeThirdDerivativeImpl>(index, othird, unweighted); + break; + case DENSE : + computeThirdDerivativeImpl>(index, othird, unweighted); + break; + case INTERCEPT : + computeThirdDerivativeImpl>(index, othird, unweighted); + break; + } + } + +#ifdef CYCLOPS_DEBUG_TIMING +#ifndef CYCLOPS_DEBUG_TIMING_LOW + auto end = bsccs::chrono::steady_clock::now(); + ///////////////////////////" + duration["compThird "] += bsccs::chrono::duration_cast(end - start).count(); +#endif +#endif +} + +template template +void ModelSpecifics::computeThirdDerivativeImpl(int index, double *othird, Weights w) { + +#ifdef CYCLOPS_DEBUG_TIMING +#ifdef CYCLOPS_DEBUG_TIMING_LOW + auto start = bsccs::chrono::steady_clock::now(); +#endif +#endif + + RealType third = static_cast(0); + RealType hessian = static_cast(0); + + if (BaseModel::cumulativeGradientAndHessian) { // Compile-time switch + + if (sparseIndices[index] == nullptr || sparseIndices[index]->size() > 0) { + + IteratorType it(sparseIndices[index].get(), N); + + RealType accNumerPid = static_cast(0); + RealType accNumerPid2 = static_cast(0); + + // find start relavent accumulator reset point + auto reset = begin(accReset); + while( *reset < it.index() ) { + ++reset; + } + + for (; it; ) { + int i = it.index(); + + if (*reset <= i) { + accNumerPid = static_cast(0.0); + accNumerPid2 = static_cast(0.0); + ++reset; + } + + const auto numerator1 = numerPid[i]; + const auto numerator2 = numerPid2[i]; + + accNumerPid += numerator1; + accNumerPid2 += numerator2; + + // Compile-time delegation + BaseModel::incrementThirdDerivative(it, + w, // Signature-only, for iterator-type specialization + &third, accNumerPid, accNumerPid2, + accDenomPid[i], hNWeight[i], 0.0, hXBeta[i], hY[i]); // When function is in-lined, compiler will only use necessary arguments + ++it; + + if (IteratorType::isSparse) { + + const int next = it ? it.index() : N; + for (++i; i < next; ++i) { + + if (*reset <= i) { + accNumerPid = static_cast(0.0); + accNumerPid2 = static_cast(0.0); + ++reset; + } + + BaseModel::incrementThirdDerivative(it, + w, // Signature-only, for iterator-type specialization + &third, accNumerPid, accNumerPid2, + accDenomPid[i], hNWeight[i], static_cast(0), hXBeta[i], hY[i]); // When function is in-lined, compiler will only use necessary arguments + } + } + } + } else { + throw new std::logic_error("Not yet support"); + } + } + + *othird = static_cast(third); + +#ifdef CYCLOPS_DEBUG_TIMING +#ifdef CYCLOPS_DEBUG_TIMING_LOW + auto end = bsccs::chrono::steady_clock::now(); + ///////////////////////////" + auto name = "compThird" + IteratorType::name + " "; + duration[name] += bsccs::chrono::duration_cast(end - start).count(); +#endif +#endif +} template void ModelSpecifics::computeFisherInformation(int indexOne, int indexTwo, diff --git a/src/cyclops/io/OutputWriter.h b/src/cyclops/io/OutputWriter.h index b9328912..89da1dff 100644 --- a/src/cyclops/io/OutputWriter.h +++ b/src/cyclops/io/OutputWriter.h @@ -234,6 +234,7 @@ class DiagnosticsOutputWriter : public BaseOutputWriter case MAX_ITERATIONS : return "MAX_ITERATIONS"; case ILLCONDITIONED : return "ILLCONDITIONED"; case MISSING_COVARIATES : return "MISSING_COVARIATES"; + case POOR_BLR_STEP : return "POOR_BLR_STEP"; default : return "FAILED"; } } diff --git a/src/cyclops/priors/CovariatePrior.cpp b/src/cyclops/priors/CovariatePrior.cpp index a1e47f1a..3ab10663 100644 --- a/src/cyclops/priors/CovariatePrior.cpp +++ b/src/cyclops/priors/CovariatePrior.cpp @@ -1,11 +1,27 @@ #include #include "priors/CovariatePrior.h" +#include "CyclicCoordinateDescent.h" namespace bsccs { namespace priors { - double FusedLaplacePrior::logDensity(const DoubleVector& betaVector, const int index) const { + + double JeffreysPrior::logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { + + double hessian = ccd.getHessianDiagonal(index); + double penalty = 0.5 * std::log(std::abs(hessian)); + return penalty; + } + + double JeffreysPrior::getDelta(GradientHessian gh,const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { + + double jerk = ccd.getJerkDiagonal(index); + //std::cerr << "j: " << jerk << std::endl; + return -((gh.first - jerk) / gh.second); + } + + double FusedLaplacePrior::logDensity(const DoubleVector& betaVector, const int index, CyclicCoordinateDescent& ccd) const { auto x = betaVector[index]; auto lambda = getLambda(); auto epsilon = getEpsilon(); @@ -41,7 +57,7 @@ namespace bsccs { } } // namespace details - double FusedLaplacePrior::getDelta(const GradientHessian gh, const DoubleVector& betaVector, const int index) const { + double FusedLaplacePrior::getDelta(const GradientHessian gh, const DoubleVector& betaVector, const int index, CyclicCoordinateDescent& ccd) const { const auto t1 = getLambda(); const auto t2 = getEpsilon(); const auto beta = betaVector[index]; @@ -222,13 +238,13 @@ namespace bsccs { // } // } - double HierarchicalNormalPrior::logDensity(const DoubleVector& beta, const int index) const { + double HierarchicalNormalPrior::logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { auto x = beta[index]; double sigma2Beta = getVariance(); return -0.5 * std::log(2.0 * PI * sigma2Beta) - 0.5 * x * x / sigma2Beta; } - double HierarchicalNormalPrior::getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index) const { + double HierarchicalNormalPrior::getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index, CyclicCoordinateDescent& ccd) const { double sigma2Beta = getVariance(); double beta = betaVector[index]; return - (gh.first + (beta / sigma2Beta)) / @@ -249,6 +265,10 @@ namespace bsccs { break; case BAR_UPDATE : prior = bsccs::make_shared(variance); + break; + case JEFFREYS : + prior = bsccs::make_shared(); + break; default : break; } return prior; diff --git a/src/cyclops/priors/CovariatePrior.h b/src/cyclops/priors/CovariatePrior.h index 36d75e32..bd17feda 100644 --- a/src/cyclops/priors/CovariatePrior.h +++ b/src/cyclops/priors/CovariatePrior.h @@ -23,6 +23,9 @@ #include "Types.h" namespace bsccs { + +class CyclicCoordinateDescent; // forward reference + namespace priors { @@ -138,9 +141,10 @@ class CovariatePrior { virtual const std::string getDescription() const = 0; // pure virtual - virtual double logDensity(const DoubleVector& beta, const int index) const = 0; // pure virtual + virtual double logDensity(const DoubleVector& beta, const int index, + CyclicCoordinateDescent& ccd) const = 0; // pure virtual - virtual double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index) const = 0; // pure virtual + virtual double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const = 0; // pure virtual virtual bool getIsRegularized() const = 0; // pure virtual @@ -175,7 +179,7 @@ class NoPrior : public CovariatePrior { return "None"; } - double logDensity(const DoubleVector& beta, const int index) const { + double logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { return 0.0; } @@ -183,7 +187,7 @@ class NoPrior : public CovariatePrior { return false; } - double getDelta(GradientHessian gh,const DoubleVector& beta, const int index) const { + double getDelta(GradientHessian gh,const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { return -(gh.first / gh.second); // No regularization } @@ -196,6 +200,46 @@ class NoPrior : public CovariatePrior { } }; +class JeffreysPrior : public CovariatePrior { +public: + JeffreysPrior() : CovariatePrior() { } + + virtual ~JeffreysPrior() { } + + const std::string getDescription() const { + return "Jeffreys (univariable)"; + } + + double logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const; + + // double logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { + // double gradient, negativeHessian; + // + // double hessian = ccd.getHessianDiagonal(index); + // + // // modelSpecifics.updateXBeta(0.0, index, false); + // // modelSpecifics.computeRemainingStatistics(false); + // // modelSpecifics.computeGradientAndHessian(index, &gradient, &negativeHessian, false); + // + // + // double penalty = -0.5 * std::log(std::abs(hessian)); + // std::cerr << hessian << " " << penalty << std::endl; + // returnā‰„ penalty; + // } + + double getDelta(GradientHessian gh,const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const; + + bool getIsRegularized() const { return true; } + + bool getSupportsKktSwindle() const { return false; } + + double getKktBoundary() const { return 0.0; } + + std::vector getVarianceParameters() const { + return std::vector(); + } +}; + class LaplacePrior : public CovariatePrior { public: @@ -218,7 +262,7 @@ class LaplacePrior : public CovariatePrior { info << "Laplace(" << lambda << ")"; return info.str(); } - double logDensity(const DoubleVector& beta, const int index) const { + double logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { auto x = beta[index]; auto lambda = getLambda(); return std::log(0.5 * lambda) - lambda * std::abs(x); @@ -237,7 +281,7 @@ class LaplacePrior : public CovariatePrior { return lambda; } - double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index) const { + double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index, CyclicCoordinateDescent& ccd) const { double beta = betaVector[index]; double lambda = getLambda(); @@ -353,9 +397,9 @@ class FusedLaplacePrior : public LaplacePrior { return info.str(); } - double logDensity(const DoubleVector& beta, const int index) const; + double logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const; - double getDelta(const GradientHessian gh, const DoubleVector& betaVector, const int index) const; + double getDelta(const GradientHessian gh, const DoubleVector& betaVector, const int index, CyclicCoordinateDescent& ccd) const; private: double getEpsilon() const { @@ -401,13 +445,13 @@ class NormalPrior : public CovariatePrior { return 0.0; } - double logDensity(const DoubleVector& beta, const int index) const { + double logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { auto x = beta[index]; double sigma2Beta = getVariance(); return -0.5 * std::log(2.0 * PI * sigma2Beta) - 0.5 * x * x / sigma2Beta; } - double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index) const { + double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index, CyclicCoordinateDescent& ccd) const { double sigma2Beta = getVariance(); double beta = betaVector[index]; return - (gh.first + (beta / sigma2Beta)) / @@ -480,13 +524,13 @@ class BarUpdatePrior : public CovariatePrior { return 0.0; } - double logDensity(const DoubleVector& beta, const int index) const { + double logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { auto x = beta[index]; double sigma2Beta = getVariance(); return -0.5 * std::log(2.0 * PI * sigma2Beta) - 0.5 * x * x / sigma2Beta; } - double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index) const { + double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index, CyclicCoordinateDescent& ccd) const { double beta = betaVector[index]; double lambda = 1 / getVariance(); @@ -558,9 +602,9 @@ class HierarchicalNormalPrior : public NormalPrior { virtual ~HierarchicalNormalPrior() { } - double logDensity(const DoubleVector& beta, const int index) const; + double logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const; - double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index) const; + double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index, CyclicCoordinateDescent& ccd) const; std::vector getVarianceParameters() const { auto tmp = NormalPrior::getVarianceParameters(); diff --git a/src/cyclops/priors/JointPrior.h b/src/cyclops/priors/JointPrior.h index 86d09141..3ad6939c 100644 --- a/src/cyclops/priors/JointPrior.h +++ b/src/cyclops/priors/JointPrior.h @@ -29,9 +29,9 @@ class JointPrior { // virtual std::vector getVariance() const = 0; // pure virtual - virtual double logDensity(const DoubleVector& beta) const = 0; // pure virtual + virtual double logDensity(const DoubleVector& beta, CyclicCoordinateDescent& ccd) const = 0; // pure virtual - virtual double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index) const = 0; // pure virtual + virtual double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const = 0; // pure virtual virtual const std::string getDescription() const = 0; // pure virtual @@ -120,12 +120,12 @@ class MixtureJointPrior : public JointPrior { // return std::move(tmp); // } - double logDensity(const DoubleVector& beta) const { + double logDensity(const DoubleVector& beta, CyclicCoordinateDescent& ccd) const { // TODO assert(beta.size() == listPriors.size()); double result = 0.0; for (size_t i = 0; i < beta.size(); ++i) { - result += listPriors[i]->logDensity(beta, i); + result += listPriors[i]->logDensity(beta, i, ccd); } return result; } @@ -134,8 +134,8 @@ class MixtureJointPrior : public JointPrior { return listPriors[index]->getIsRegularized(); } - double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index) const { - return listPriors[index]->getDelta(gh, beta, index); + double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { + return listPriors[index]->getDelta(gh, beta, index, ccd); } bool getSupportsKktSwindle(const int index) const { @@ -506,11 +506,11 @@ class HierarchicalJointPrior : public JointPrior { // return hierarchyPriors[level]->getVariance(0); // } - double logDensity(const DoubleVector& beta) const { + double logDensity(const DoubleVector& beta, CyclicCoordinateDescent& ccd) const { double result = 0.0; // for (DoubleVector::const_iterator it = beta.begin(); it != beta.end(); ++it) { for (size_t i = 0; i < beta.size(); ++i) { - result += hierarchyPriors[0]->logDensity(beta, i); + result += hierarchyPriors[0]->logDensity(beta, i, ccd); } return result; } @@ -531,7 +531,7 @@ class HierarchicalJointPrior : public JointPrior { return 0.0; // TODO fix } - double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index) const { + double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { // double t1 = 1/hierarchyPriors[0]->getVariance(0); // this is the hyperparameter that is used in the original code // double t2 = 1/hierarchyPriors[1]->getVariance(0); @@ -604,17 +604,17 @@ class FullyExchangeableJointPrior : public JointPrior { // return variances; // } - double logDensity(const DoubleVector& beta) const { + double logDensity(const DoubleVector& beta, CyclicCoordinateDescent& ccd) const { double result = 0.0; // for (DoubleVector::const_iterator it = beta.begin(); it != beta.end(); ++it) { for (size_t i = 0; i < beta.size(); ++i) { - result += singlePrior->logDensity(beta, i); + result += singlePrior->logDensity(beta, i, ccd); } return result; } - double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index) const { - return singlePrior->getDelta(gh, beta, index); + double getDelta(const GradientHessian gh, const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { + return singlePrior->getDelta(gh, beta, index, ccd); } bool getIsRegularized(const int index) const { diff --git a/src/cyclops/priors/NewCovariatePrior.h b/src/cyclops/priors/NewCovariatePrior.h index c27f13bf..6999a9b0 100644 --- a/src/cyclops/priors/NewCovariatePrior.h +++ b/src/cyclops/priors/NewCovariatePrior.h @@ -47,7 +47,7 @@ class NewLaplacePrior : public CovariatePrior { info << "Laplace(" << locationLambda.first << ", " << locationLambda.second << ")"; return info.str(); } - double logDensity(const DoubleVector& beta, const int index) const { + double logDensity(const DoubleVector& beta, const int index, CyclicCoordinateDescent& ccd) const { const auto x = beta[index]; const auto locationLambda = getLocationLambda(); const auto location = locationLambda.first; @@ -69,7 +69,7 @@ class NewLaplacePrior : public CovariatePrior { return lambda; } - double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index) const { + double getDelta(GradientHessian gh, const DoubleVector& betaVector, const int index, CyclicCoordinateDescent& ccd) const { const auto locationLambda = getLocationLambda(); const double location = locationLambda.first; @@ -184,6 +184,9 @@ static PriorPtr makePrior(PriorType priorType, PriorFunctionPtr& priorFunction, case BAR_UPDATE : Rcpp::stop("Parameterized BAR updates are not yet implemented"); break; + case JEFFREYS : + Rcpp::stop("Parameterized Jeffreys priors are not yet implemented"); + break; default : break; } return prior; diff --git a/tests/testthat/test-RcppDataModel.R b/tests/testthat/test-RcppDataModel.R new file mode 100644 index 00000000..041a0c10 --- /dev/null +++ b/tests/testthat/test-RcppDataModel.R @@ -0,0 +1,41 @@ +library(testthat) + +test_that("quantile", { + set.seed(123) + vector <- runif(100) + + expect_equal( + Cyclops:::.cyclopsQuantile(vector, q = 0.5), + median(vector)) + expect_equal( + Cyclops:::.cyclopsMedian(vector), + median(vector)) + + expect_error( + Cyclops:::.cyclopsQuantile(vector, q = -0.1), + "Invalid quantile") +}) + +test_that("null ptr", { + expect_error( + Cyclops:::.isRcppPtrNull(1), + "Input must be an Rcpp externalptr") +}) + +test_that("print MatrixMarket format", { + binomial_bid <- c(1,5,10,20,30,40,50,75,100,150,200) + binomial_n <- c(31,29,27,25,23,21,19,17,15,15,15) + binomial_y <- c(0,3,6,7,9,13,17,12,11,14,13) + + log_bid <- log(c(rep(rep(binomial_bid, binomial_n - binomial_y)), rep(binomial_bid, binomial_y))) + y <- c(rep(0, sum(binomial_n - binomial_y)), rep(1, sum(binomial_y))) + data <- createCyclopsData(y ~ log_bid, modelType = "lr") + + tempFile <- tempfile() + Cyclops:::printMatrixMarket(data, file = tempFile) + text <- readr::read_delim(file = tempFile, delim = " ", skip = 3, col_names = FALSE, col_types = readr::cols()) + expect_equal(nrow(text), getNumberOfCovariates(data) * getNumberOfRows(data)) + expect_equal(ncol(text), 1 + getNumberOfCovariates(data)) + + unlink(tempFile) +}) diff --git a/tests/testthat/test-covariateRegularization.R b/tests/testthat/test-covariateRegularization.R index 26ada2d1..74960085 100644 --- a/tests/testthat/test-covariateRegularization.R +++ b/tests/testthat/test-covariateRegularization.R @@ -104,7 +104,7 @@ test_that("Preclude intercept regularization by default", { expect_equal(coef(c2), coef(c3)) - expect_less_than(coef(c1)[1], # Intercept is regularized + expect_lt(coef(c1)[1], # Intercept is regularized coef(c2)[1]) }) @@ -124,6 +124,14 @@ test_that("Specify each prior independently", { expect_true(coef(cyclopsFit)[4] == 0) expect_true(coef(cyclopsFit)[5] != 0) + + expect_equal(getHyperParameter(cyclopsFit), c(1,1)) + + expect_false(Cyclops:::.cyclopsGetIsRegularized(cyclopsFit$interface, 0)) + expect_true(Cyclops:::.cyclopsGetIsRegularized(cyclopsFit$interface, 1)) + + expect_equal(Cyclops:::.cyclopsGetLogLikelihood(cyclopsFit$interface), + cyclopsFit$log_likelihood) }) diff --git a/tests/testthat/test-cv.R b/tests/testthat/test-cv.R index 1ac28654..8b47229c 100644 --- a/tests/testthat/test-cv.R +++ b/tests/testthat/test-cv.R @@ -199,10 +199,10 @@ test_that("Seed gets returned", { y <- 0 x <- 1 data <- createCyclopsData(y ~ x, modelType = "lr") - fit <- fitCyclopsModel(data, control = createControl(seed = 123)) + fit <- fitCyclopsModel(data, control = createControl(seed = 123), warnings = FALSE) expect_equal(fit$seed, 123) fit <- fitCyclopsModel(data, forceNewObject = TRUE, - control = createControl(seed = NULL)) + control = createControl(seed = NULL), warnings = FALSE) expect_true(!is.null(fit$seed)) }) diff --git a/tests/testthat/test-fastBarUpdate.R b/tests/testthat/test-fastBarUpdate.R index d54b76f6..482bed9c 100644 --- a/tests/testthat/test-fastBarUpdate.R +++ b/tests/testthat/test-fastBarUpdate.R @@ -41,7 +41,7 @@ test_that("Small Poisson dense regression; unpenalized regression", { prior <- createPrior(priorType = "barupdate", variance = Inf) control <- createControl(convergenceType = "onestep") fit <- fitCyclopsModel(data, prior = prior, control = control, - startingCoefficients = b.old) + startingCoefficients = b.old, warnings = FALSE) delta <- max(abs((coef(fit) - b.old))) b.old <- coef(fit) } diff --git a/tests/testthat/test-finiteMLE.R b/tests/testthat/test-finiteMLE.R new file mode 100644 index 00000000..ab76f8d5 --- /dev/null +++ b/tests/testthat/test-finiteMLE.R @@ -0,0 +1,103 @@ +library("testthat") +library("survival") + +context("test-infiniteMLE.R") + +simulate <- function() { + set.seed(123) + n <- 1000 + timeToCensor <- rexp(n, 0.01) + exposure <- runif(n) < 0.5 + timeToEvent <- rep(Inf, n) + timeToEvent[!exposure] <- rexp(sum(!exposure), 0.01) + time <- pmin(timeToCensor, timeToEvent) + outcome <- timeToEvent < timeToCensor + list(time = time, outcome = outcome, exposure = exposure) +} + +test_that("Check for infinite MLE in Cox example with no outcomes in one treatment arm", { + + data <- simulate() + + expect_warning(gold <- coxph(Surv(time, outcome) ~ exposure, data = data), + regexp = ".*coefficient may be infinite.*") + + 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_true(is.na(ci[2])) +}) + +test_that("Check multivariable Jeffreys prior", { + + data <- simulate() + data$x2 <- runif(length(data$time), 0.5) < 0.5 + + cyclopsData <- createCyclopsData(Surv(time, outcome) ~ exposure + x2, data = data, modelType = "cox") + expect_error(fit <- fitCyclopsModel(cyclopsData, + prior = createPrior(priorType = "jeffreys")), + regexp = ".*1 covariate.*") +}) + +test_that("Check Jeffreys prior with correct control", { + + data <- simulate() + + cyclopsData <- createCyclopsData(Surv(time, outcome) ~ exposure, data = data, modelType = "cox") + expect_warning(fit <- fitCyclopsModel(cyclopsData, + prior = createPrior(priorType = "jeffreys")), + regexp = "BLR convergence") + + fit <- fitCyclopsModel(cyclopsData = cyclopsData, forceNewObject = TRUE, + prior = createPrior(priorType = "jeffreys"), + control = createControl(convergenceType = "lange")) + expect_true(fit$return_flag == "SUCCESS") +}) + +test_that("Check Jeffreys prior with indicator only", { + + data <- data.frame(time = c(1,2), + outcome = c(1,0), + exposure = c(1,2)) + + cyclopsData <- createCyclopsData(Surv(time, outcome) ~ exposure, data = data, modelType = "cox") + expect_error(fitCyclopsModel(cyclopsData, prior = createPrior(priorType = "jeffreys")), + regexp = "*.for indicator covariates") + +}) + +test_that("Check univariable Jeffreys prior", { + + data <- simulate() + + cyclops1 <- createCyclopsData(Surv(time, outcome) ~ exposure, data = data, modelType = "cox") + cyclops2 <- createCyclopsData(Surv(time, outcome) ~ exposure, data = data, modelType = "cox") + + expect_warning(uniform <- fitCyclopsModel(cyclops1)) + + jeffreys <- fitCyclopsModel(cyclops2, + prior = createPrior(priorType = "jeffreys"), + control = createControl(convergenceType = "lange", + tolerance = 1E-7), + forceNewObject = TRUE) + + jci <- confint(jeffreys, parm = "exposureTRUE", level = 0.95, overrideNoRegularization = TRUE) + expect_false(is.na(jci[2])) + + # grid <- seq(from=-8, to = -2, length.out = 1000) + # uniformProfile <- getCyclopsProfileLogLikelihood(uniform, parm = "exposureTRUE", x = grid) + # jeffreysProfile <- getCyclopsProfileLogLikelihood(jeffreys, parm = "exposureTRUE", x = grid) + # + # plot(uniformProfile$point, uniformProfile$value - max(uniformProfile$value), type = "l") + # lines(jeffreysProfile$point, jeffreysProfile$value - max(uniformProfile$value, lwd = 2)) + # + # plot(jeffreysProfile$point, jeffreysProfile$value - max(jeffreysProfile$value), type = "l") + # + # maxProfile <- function(profile) { + # idx <- which(profile$value == max(profile$value)) + # list(point = profile$point[idx], value = profile$value[idx]) + # } + # + # maxProfile(jeffreysProfile) +}) diff --git a/tests/testthat/test-int64.R b/tests/testthat/test-int64.R index c38314ba..81b9d3eb 100644 --- a/tests/testthat/test-int64.R +++ b/tests/testthat/test-int64.R @@ -25,6 +25,18 @@ test_that("Check int64 handling", { int64Result <- coef(fit64) expect_equal(int64Result, int32Result) + + # Change to int64 labels + longName <- "1234567890" + data$covariates$covariateId[which(data$covariates$covariateId == bit64::as.integer64(1))] <- bit64::as.integer64(longName) + cyclopsData64b <- convertToCyclopsData(data$outcomes,data$covariates, modelType = "lr", addIntercept = TRUE) + fit2 <- fitCyclopsModel(cyclopsData64b) + result2 <- coef(fit2) + expect_error( + expect_equal(sort(result2), sort(int64Result)), + "1 string mismatch") + expect_equal(rownames(confint(fit2, parm = longName)), longName) + expect_equal(rownames(confint(fit2, parm = bit64::as.integer64(longName))), longName) }) test_that("Check int64 reductions", { diff --git a/tests/testthat/test-multitypePoisson.R b/tests/testthat/test-multitypePoisson.R index 6ea0b267..e565b493 100644 --- a/tests/testthat/test-multitypePoisson.R +++ b/tests/testthat/test-multitypePoisson.R @@ -18,8 +18,8 @@ test_that("Small multi-type Poisson dense regression", { goldFit <- glm(counts ~ outcome + treatment, data = dobson1, family = poisson()) - glmFit <- glm(counts ~ outcome + treatment, data = dobson, contrasts = dobson$type, - family = poisson()) # gold standard + # glmFit <- glm(counts ~ outcome + treatment, data = dobson, contrasts = dobson$type, + # family = poisson()) # gold standard dataPtrD <- createCyclopsData(counts ~ outcome + treatment, data = dobson, type = dobson$type, @@ -125,8 +125,8 @@ test_that("Small multi-type Poisson with hierarchical prior", { dobson$type = as.factor(c(rep("A",9),rep("B",9))) tolerance <- 1E-4 - glmFit <- glm(counts ~ outcome + treatment, data = dobson, contrasts = dobson$type, - family = poisson()) # gold standard + # glmFit <- glm(counts ~ outcome + treatment, data = dobson, contrasts = dobson$type, + # family = poisson()) # gold standard dataPtrD <- createCyclopsData(Multitype(counts, type) ~ outcome + treatment, data = dobson, modelType = "pr") diff --git a/tests/testthat/test-parameterizedPriors.R b/tests/testthat/test-parameterizedPriors.R index e8c46b56..9194311c 100644 --- a/tests/testthat/test-parameterizedPriors.R +++ b/tests/testthat/test-parameterizedPriors.R @@ -66,7 +66,8 @@ test_that("Specify parameterized L1 regularization", { fit <- fitCyclopsModel(data, prior, forceNewObject = TRUE), "Excluding intercept") - comp <- fitCyclopsModel(data, prior = createPrior("laplace", variance = 2), forceNewObject = TRUE) + comp <- fitCyclopsModel(data, prior = createPrior("laplace", variance = 2), + warnings = FALSE, forceNewObject = TRUE) expect_equal(coef(fit), coef(comp)) @@ -78,7 +79,7 @@ test_that("Specify parameterized L1 regularization", { }) }) - fit <- fitCyclopsModel(data, prior, forceNewObject = TRUE) + fit <- fitCyclopsModel(data, prior, warnings = FALSE, forceNewObject = TRUE) expect_equivalent(coef(fit)[4:5], c(0.2, 0.2)) }) @@ -126,7 +127,7 @@ test_that("Using parameterized cross-validation", { if (i < (ncovars / 2)) { return(c(0,x[1])) } else { - return(c(0,2 * x[2])) + return(c(0,x[2])) } } ) @@ -134,11 +135,25 @@ test_that("Using parameterized cross-validation", { values = c(1, 0.5), useCrossValidation = TRUE) - # time3 <- system.time(fit3 <- fitCyclopsModel(cyclopsData, - # prior=prior3, - # control=control, - # forceNewObject = TRUE)) + control3 <- createAutoGridCrossValidationControl(outerGrid = c(0.1, 0.1, 0.2)) + time3 <- system.time(fit3 <- fitCyclopsModel(cyclopsData, + prior=prior3, + control=control3, + forceNewObject = TRUE)) + + expect_equal(length(unlist(strsplit(fit3$cross_validation, split = " "))), 3) + + expect_error(createAutoGridCrossValidationControl(outerGrid = c(1,2,3), + autoPosition = 0), + "Auto-position is invalid") + # The following should execute w/o error + createAutoGridCrossValidationControl(outerGrid = matrix(c(1,2,3,4,5,6), nrow = 2), + autoPosition = 2) + createAutoGridCrossValidationControl(outerGrid = matrix(c(1,2,3,4,5,6), nrow = 2), + autoPosition = 3) }) + + diff --git a/tests/testthat/test-profileLikelihood.R b/tests/testthat/test-profileLikelihood.R index 9918a3ca..1716d4d1 100644 --- a/tests/testthat/test-profileLikelihood.R +++ b/tests/testthat/test-profileLikelihood.R @@ -25,6 +25,28 @@ start, length, event, x1, x2 expect_equivalent(coef(fit)["x1"], argMax, tolerance = 0.01) }) +test_that("Check evaluate profile likelihood with multiple threads", { + test <- read.table(header=T, sep = ",", text = " +start, length, event, x1, x2 +0, 4, 1,0,0 +0, 3.5,1,2,0 +0, 3, 0,0,1 +0, 2.5,1,0,1 +0, 2, 1,1,1 +0, 1.5,0,1,0 +0, 1, 1,1,0 +") + data <- createCyclopsData(Surv(length, event) ~ x1 + x2, data = test, + modelType = "cox") + fit <- fitCyclopsModel(data, control = createControl(threads = 2)) + + x <- seq(from = 0.8, to = 0.9, length = 100) + out <- getCyclopsProfileLogLikelihood(fit, "x1", x) + + argMax <- out$point[which(out$value == max(out$value))] + expect_equivalent(coef(fit)["x1"], argMax, tolerance = 0.01) +}) + test_that("Check evaluate profile likelihood with one variable; cold-start", { test <- read.table(header=T, sep = ",", text = " start, length, event, x1, x2 diff --git a/tests/testthat/test-simulation.R b/tests/testthat/test-simulation.R new file mode 100644 index 00000000..db9cb7a1 --- /dev/null +++ b/tests/testthat/test-simulation.R @@ -0,0 +1,37 @@ +library("testthat") + +context("test-simulation.R") + +# +# Simulation +# + +test_that("Check logistic output from other programs", { + set.seed(123) + tolerance <- 1E-2 + data <- simulateCyclopsData(nstrata = 1, nrows = 1000, ncovars = 2, zeroEffectSizeProp = 0.0, eCovarsPerRow = 0.5, + model = "logistic") + suppressWarnings(cyclopsFit <- fitCyclopsSimulation(data, model = "logistic", useCyclops = TRUE)) + otherFit <- fitCyclopsSimulation(data, model = "logistic", useCyclops = FALSE) + expect_equal(cyclopsFit$coef, otherFit$coef, tolerance = tolerance) +}) + +test_that("Check survival output from other programs", { + set.seed(123) + tolerance <- 1E-2 + data <- simulateCyclopsData(nstrata = 1, nrows = 1000, ncovars = 2, zeroEffectSizeProp = 0.0, eCovarsPerRow = 0.5, + model = "survival") + suppressWarnings(cyclopsFit <- fitCyclopsSimulation(data, model = "survival", useCyclops = TRUE)) + otherFit <- fitCyclopsSimulation(data, model = "survival", useCyclops = FALSE) + expect_equal(cyclopsFit$coef, otherFit$coef, tolerance = tolerance) + + plot <- Cyclops:::plotCyclopsSimulationFit(cyclopsFit, log(data$effectSizes$rr), label = "Cyclops") + expect_false(is.null(plot)) +}) + +test_that("Check simulation edge cases", { + expect_error(simulateCyclopsData(model = "not_known"), + "Unknown model") + + expect_equal(mse(goldStandard = c(0,2), estimates = c(0,0)), 2) +}) diff --git a/tests/testthat/test-smallBernoulli.R b/tests/testthat/test-smallBernoulli.R index efd3e5d1..08380993 100644 --- a/tests/testthat/test-smallBernoulli.R +++ b/tests/testthat/test-smallBernoulli.R @@ -25,6 +25,32 @@ test_that("Small Bernoulli dense regression", { expect_equal(cyclopsFitD$log_likelihood, logLik(glmFit)[[1]], tolerance = tolerance) expect_equal(confint(cyclopsFitD, c(1:2))[,2:3], confint(glmFit, c(1:2)), tolerance = tolerance) expect_equal(predict(cyclopsFitD), predict(glmFit, type = "response"), tolerance = tolerance) + + # Try with Zhang-Oles criterion + zoFit <- fitCyclopsModel(dataPtrD, prior = createPrior("none"), forceNewObject = TRUE, + control = createControl(convergenceType = "zhang", noiseLevel = "noisy")) + expect_equal(coef(cyclopsFitD), coef(zoFit), tolerance = tolerance) + + # Try with KKT swindle + gold <- fitCyclopsModel(dataPtrD, prior = createPrior("laplace"), forceNewObject = TRUE, warnings = FALSE) + + kktFit <- fitCyclopsModel(dataPtrD, prior = createPrior("laplace"), forceNewObject = TRUE, warnings = FALSE, + control = createControl(useKKTSwindle = TRUE, noiseLevel = "quiet")) + expect_equal(coef(gold), coef(kktFit)) + + # Test old logging routines + file <- tempfile() + Cyclops:::.cyclopsLogResults(cyclopsFitD$interface, file, TRUE) + result <- read.csv(file) + expect_equivalent(result$estimate, coef(cyclopsFitD), tolerance = tolerance) + + cap <- capture.output(print(dataPtrD)) + expect_false(is.null(cap)) +}) + +test_that("Old BLR format", { + data <- readCyclopsData(system.file("extdata/infert_ccd.txt", package = "Cyclops"), "clr") + expect_equal(getNumberOfRows(data), 248) }) test_that("Add intercept via finalize", { diff --git a/tests/testthat/test-smallCox.R b/tests/testthat/test-smallCox.R index 9c4f6e20..873c5b77 100644 --- a/tests/testthat/test-smallCox.R +++ b/tests/testthat/test-smallCox.R @@ -182,10 +182,18 @@ start, length, event, x1, x2 weights = rep(c(0,1), 7)) weights <- rep(c(1,0), 7) - # predictiveLogLik <- getCyclopsPredictiveLogLikelihood(cyclopsFitClean, weights) - # - # expect_equal(predictiveLogLik, logLik(goldRight)[[1]]) + predictiveLogLik <- Cyclops:::getCyclopsPredictiveLogLikelihood(cyclopsFitClean, weights) + expect_equal(predictiveLogLik, logLik(goldRight)[[1]]) + # Test error cases in predictive likelihood + expect_error(Cyclops:::getCyclopsPredictiveLogLikelihood(cyclopsFitClean, weights = c(1,2)), + "Must provide") + + expect_error(Cyclops:::getCyclopsPredictiveLogLikelihood(cyclopsFitClean, weights = rep(-1, 14)), + "Only non-negative weights") + + cap <- capture.output(print(cyclopsFitClean)) + expect_false(is.null(cap)) }) test_that("Check SQL interface for a very small Cox example with failure ties and strata", { diff --git a/tests/testthat/test-smallFineGray.R b/tests/testthat/test-smallFineGray.R index 2ceeb0bb..5e2316f8 100644 --- a/tests/testthat/test-smallFineGray.R +++ b/tests/testthat/test-smallFineGray.R @@ -252,7 +252,7 @@ test_that("Check very small Fine-Gray versus Cox for 0/1 events (censorWeights d dataPtr <- Cyclops:::createCyclopsData(fgDat$surv ~ test$x1 + test$x2, modelType = "fgr", censorWeights = fgDat$weights) cyclopsFit1 <- Cyclops:::fitCyclopsModel(dataPtr) dataPtr <- Cyclops:::createCyclopsData(fgDat$surv ~ test$x1 + test$x2, modelType = "cox", censorWeights = fgDat$weights) - cyclopsFit2 <- Cyclops:::fitCyclopsModel(dataPtr) + cyclopsFit2 <- Cyclops:::fitCyclopsModel(dataPtr, warnings = FALSE) tolerance <- 1E-4 expect_equivalent(coef(cyclopsFit1), coef(cyclopsFit2), tolerance = tolerance) @@ -272,7 +272,7 @@ test_that("Check very small Cox example with weights and censor weights defined fgDat <- Cyclops:::getFineGrayWeights(test$length, test$event) dataPtr <- Cyclops:::createCyclopsData(fgDat$surv ~ test$x1 + test$x2, modelType = "cox", weights = c(1, 1, 1, 0, 0, 1, 1), censorWeights = fgDat$weights) - cyclopsFit1 <- Cyclops:::fitCyclopsModel(dataPtr) + cyclopsFit1 <- Cyclops:::fitCyclopsModel(dataPtr, warnings = FALSE) dataPtr <- Cyclops:::createCyclopsData(fgDat$surv ~ test$x1 + test$x2, modelType = "cox", weights = c(1, 1, 1, 0, 0, 1, 1)) cyclopsFit2 <- Cyclops:::fitCyclopsModel(dataPtr) diff --git a/tests/testthat/test-smallMM.R b/tests/testthat/test-smallMM.R new file mode 100644 index 00000000..178b4708 --- /dev/null +++ b/tests/testthat/test-smallMM.R @@ -0,0 +1,62 @@ +library("testthat") + +context("test-smallMM.R") + +# +# Small MM regression +# + + +test_that("Small Bernoulli dense regression with MM algorithm", { + binomial_bid <- c(1,5,10,20,30,40,50,75,100,150,200) + binomial_n <- c(31,29,27,25,23,21,19,17,15,15,15) + binomial_y <- c(0,3,6,7,9,13,17,12,11,14,13) + + log_bid <- log(c(rep(rep(binomial_bid, binomial_n - binomial_y)), rep(binomial_bid, binomial_y))) + y <- c(rep(0, sum(binomial_n - binomial_y)), rep(1, sum(binomial_y))) + + tolerance <- 1E-3 + + gold <- glm(y ~ log_bid, family = binomial()) # gold standard + + data <- createCyclopsData(y ~ log_bid, modelType = "lr") + fit <- fitCyclopsModel(data, prior = createPrior("none"), + control = createControl(algorithm = "mm")) + expect_equal(coef(fit), coef(gold), tolerance = tolerance) +}) + +test_that("Small Bernoulli sparse regression with MM algorithm", { + binomial_bid <- c(1,5,10,20,30,40,50,75,100,150,200) + binomial_n <- c(31,29,27,25,23,21,19,17,15,15,15) + binomial_y <- c(0,3,6,7,9,13,17,12,11,14,13) + + log_bid <- log(c(rep(rep(binomial_bid, binomial_n - binomial_y)), rep(binomial_bid, binomial_y))) + y <- c(rep(0, sum(binomial_n - binomial_y)), rep(1, sum(binomial_y))) + + tolerance <- 1E-3 + + gold <- glm(y ~ log_bid, family = binomial()) # gold standard + + data <- createCyclopsData(y ~ 1, sparseFormula = ~ log_bid, modelType = "lr") + fit <- fitCyclopsModel(data, prior = createPrior("none"), + control = createControl(algorithm = "mm")) + expect_equal(coef(fit), coef(gold), tolerance = tolerance) +}) + +test_that("Small Bernoulli indicator regression with MM algorithm", { + binomial_bid <- c(0,0,0,0,0,0,1,1,1,1,1) + binomial_n <- c(31,29,27,25,23,21,19,17,15,15,15) + binomial_y <- c(0,3,6,7,9,13,17,12,11,14,13) + + log_bid <- c(rep(rep(binomial_bid, binomial_n - binomial_y)), rep(binomial_bid, binomial_y)) + y <- c(rep(0, sum(binomial_n - binomial_y)), rep(1, sum(binomial_y))) + + tolerance <- 1E-3 + + gold <- glm(y ~ log_bid, family = binomial()) # gold standard + + data <- createCyclopsData(y ~ 1, indicatorFormula = ~ log_bid, modelType = "lr") + fit <- fitCyclopsModel(data, prior = createPrior("none"), + control = createControl(algorithm = "mm")) + expect_equal(coef(fit), coef(gold), tolerance = tolerance) +}) diff --git a/tests/testthat/test-smallNormal.R b/tests/testthat/test-smallNormal.R new file mode 100644 index 00000000..e2f3d2b5 --- /dev/null +++ b/tests/testthat/test-smallNormal.R @@ -0,0 +1,37 @@ +library("testthat") + +context("test-smallNormal.R") + +# +# Small Normal regression +# + +test_that("Small Normal dense regression with CCD algorithm", { + + x <- log(c(1,5,10,20,30,40,50,75,100,150,200)) + y <- c(0,3,6,7,9,13,17,12,11,14,13) + + tolerance <- 1E-4 + + gold <- lm(y ~ x) # gold standard + + data <- createCyclopsData(y ~ x, modelType = "ls") + fit <- fitCyclopsModel(data, prior = createPrior("none"), + control = createControl(tolerance = 1E-8)) + expect_equal(coef(fit), coef(gold), tolerance = tolerance) +}) + +test_that("Small Normal sparse regression with CCD algorithm", { + + x <- log(c(1,5,10,20,30,40,50,75,100,150,200)) + y <- c(0,3,6,7,9,13,17,12,11,14,13) + + tolerance <- 1E-4 + + gold <- lm(y ~ x) # gold standard + + data <- createCyclopsData(y ~ 1, sparseFormula = ~ x, modelType = "ls") + fit <- fitCyclopsModel(data, prior = createPrior("none"), + control = createControl(tolerance = 1E-8)) + expect_equal(coef(fit), coef(gold), tolerance = tolerance) +}) diff --git a/tests/testthat/test-smallPoisson.R b/tests/testthat/test-smallPoisson.R index 6af8989c..64280568 100644 --- a/tests/testthat/test-smallPoisson.R +++ b/tests/testthat/test-smallPoisson.R @@ -26,6 +26,75 @@ test_that("Small Poisson dense regression", { expect_equal(confint(cyclopsFitD, c("(Intercept)","outcome3")), confint(cyclopsFitD, c(1,3))) }) +test_that("Small Poisson dense regression in 32bit", { + dobson <- data.frame( + counts = c(18,17,15,20,10,20,25,13,12), + outcome = gl(3,1,9), + treatment = gl(3,3) + ) + tolerance <- 1E-4 + + gold <- glm(counts ~ outcome + treatment, data = dobson, family = poisson()) # gold standard + + data <- createCyclopsData(counts ~ outcome + treatment, data = dobson, + modelType = "pr", floatingPoint = 32) + expect_equal(getFloatingPointSize(data), 32) + + fit <- fitCyclopsModel(data, + prior = createPrior("none"), + control = createControl(noiseLevel = "silent")) + + expect_equal(coef(fit), coef(gold), tolerance = tolerance) + + expect_equal(Cyclops:::.cyclopsGetMeanOffset(data), 0) +}) + +test_that("Errors for different criteria", { + dobson <- data.frame( + counts = c(18,17,15,20,10,20,25,13,12), + outcome = gl(3,1,9), + treatment = gl(3,3) + ) + tolerance <- 1E-3 + + gold <- glm(counts ~ outcome + treatment, data = dobson, family = poisson()) # gold standard + + data <- createCyclopsData(counts ~ outcome + treatment, data = dobson, + modelType = "pr") + + fit <- fitCyclopsModel(data, + prior = createPrior("none"), + control = createControl(convergenceType = "mittal", + noiseLevel = "silent")) + + expect_equal(coef(fit), coef(gold), tolerance = tolerance) +}) + +test_that("Small Poisson dense regression with offset", { + dobson <- data.frame( + counts = c(18,17,15,20,10,20,25,13,12), + outcome = as.numeric(gl(3,1,9)), + treatment = gl(3,3) + ) + tolerance <- 1E-4 + + gold <- glm(counts ~ treatment, offset = outcome, data = dobson, family = poisson()) # gold standard + + data <- createCyclopsData(counts ~ treatment, offset = as.numeric(outcome), + data = dobson, modelType = "pr") + + fit <- fitCyclopsModel(data, + prior = createPrior("none"), + control = createControl(noiseLevel = "silent")) + + expect_equal(coef(fit), coef(gold), tolerance = tolerance) + + expect_error(Cyclops:::.cyclopsGetMeanOffset(fit), + "Input must be a cyclopsData object") + + expect_equal(Cyclops:::.cyclopsGetMeanOffset(data), 2) +}) + test_that("Small Poisson fixed beta", { dobson <- data.frame( counts = c(18,17,15,20,10,20,25,13,12), diff --git a/tests/testthat/test-sorting.R b/tests/testthat/test-sorting.R new file mode 100644 index 00000000..e73d385e --- /dev/null +++ b/tests/testthat/test-sorting.R @@ -0,0 +1,16 @@ +library("testthat") + +context("test-sorting.R") + +# +# Sort functions +# + +test_that("Check sorting functions", { + set.seed(123) + vector <- c(1:100) + + expect_true(Cyclops:::.isSortedVectorList(vectorList = list(vector[1]), ascending = FALSE)) + expect_true(Cyclops:::.isSortedVectorList(vectorList = list(vector), ascending = c(TRUE))) + expect_false(Cyclops:::.isSortedVectorList(vectorList = list(vector), ascending = c(FALSE))) +}) diff --git a/tests/testthat/test-weighted_cv.R b/tests/testthat/test-weighted_cv.R index 4886458f..1d60cb4d 100644 --- a/tests/testthat/test-weighted_cv.R +++ b/tests/testthat/test-weighted_cv.R @@ -13,7 +13,7 @@ test_that("Logistic regression with cross-validation", { model="logistic") cyclopsData <- convertToCyclopsData(outcomes = sim$outcomes, covariates = sim$covariates, modelType = "lr") - fit <- fitCyclopsModel(cyclopsData = cyclopsData, + fit <- fitCyclopsModel(cyclopsData = cyclopsData, warnings = FALSE, prior = createPrior("laplace", useCrossValidation = TRUE), control = createControl(fold = 10, cvRepetitions = 1, startingVariance = 0.1, noiseLevel = "quiet", seed = 123)) diff --git a/tests/testthat/test-weighting.R b/tests/testthat/test-weighting.R index d5b3986e..73bbc64e 100644 --- a/tests/testthat/test-weighting.R +++ b/tests/testthat/test-weighting.R @@ -252,7 +252,7 @@ test_that("Small Poisson dense regression with weighting", { fit <- fitCyclopsModel(dataPtrD, prior = createPrior("laplace", useCrossValidation = TRUE), - weights = weights, + weights = weights, warnings = FALSE, control = createControl(minCVData = 1, noiseLevel = "quiet")) }) @@ -264,9 +264,11 @@ test_that("Check conditional Poisson with cross-validation",{ effectSizeSd=0.5, eCovarsPerRow=2, model="poisson") - cyclopsData <- convertToCyclopsData(outcomes = sim$outcomes, - covariates = sim$covariates, - modelType = "cpr") + suppressWarnings( + cyclopsData <- convertToCyclopsData(outcomes = sim$outcomes, + covariates = sim$covariates, + modelType = "cpr") + ) fit <- fitCyclopsModel(cyclopsData = cyclopsData, prior = createPrior("laplace", useCrossValidation = TRUE,