diff --git a/DESCRIPTION b/DESCRIPTION index 0a82571..50448ae 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: accept Title: The Acute COPD Exacerbation Prediction Tool (ACCEPT) -Version: 1.1.1 +Version: 1.2.0 Authors@R: c( person("Amin", "Adibi", email = "adibi@alumni.ubc.ca", role = c("aut", "cre")), person("Mohsen", "Sadatsafavi", email = "mohsen.sadatsafavi@ubc.ca", role = c("aut", "cph")), diff --git a/NAMESPACE b/NAMESPACE index 2a45c64..49246ff 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,6 +9,7 @@ export(accept) export(accept1) export(accept2) export(accept3) +export(accept3_cprd) export(plotExacerbations) export(plotHeatMap) export(predictCountProb) @@ -16,6 +17,10 @@ export(set_openai_api_key) export(show_openai_api_key) import(dplyr) import(tidyselect) +importFrom(dplyr,mutate) +importFrom(dplyr,select) +importFrom(dplyr,starts_with) +importFrom(dplyr,tibble) importFrom(hardhat,scream) importFrom(reldist,wtd.quantile) importFrom(splines,bs) diff --git a/R/predict.R b/R/predict.R index b2888dd..a472a41 100644 --- a/R/predict.R +++ b/R/predict.R @@ -667,7 +667,7 @@ accept2 <- function (patientData, random_sampling_N = 1e2, lastYrExacCol = "Last #' @param version indicates which version of ACCEPT needs to be called. Options include "accept1", "accept2", "accept3" (default). #' @param prediction_interval default is FALSE. If set to TRUE, returns prediction intervals of the predictions. #' @param return_predictors default is FALSE. IF set to TRUE, returns the predictors along with prediction results. -#' @param country Required for accept3 version. Three-letter ISO country code (e.g., "CAN", "USA", "GBR"). Supported countries: ARG, AUS, BRA, CAN, COL, DEU, DNK, ESP, FRA, GBR, ITA, JPN, KOR, MEX, NLD, NOR, SWE, USA. For unsupported countries, provide obs_modsev_risk parameter or add it as a column in the data. +#' @param country Required for accept3 version. Three-letter ISO country code (e.g., "CAN", "USA", "GBR"). Supported countries: ARG, AUS, BRA, CAN, COL, DEU, DNK, ESP, FRA, GBR, ITA, JPN, KOR, MEX, NLD, NOR, SWE, USA. For the UK, two care-setting-specific versions are available: "GBR-primary" (primary care, ACCEPT 3.0-CPRD) and "GBR-specialty" (specialty care). Bare "GBR" defaults to specialty care with a warning. For unsupported countries, provide obs_modsev_risk parameter or add it as a column in the data. #' @param obs_modsev_risk Observed moderate-to-severe exacerbation risk for unsupported countries. Can be provided as a parameter or as a column in newdata. #' @param ... for other versions of accept. #' @return patientData with prediction. @@ -700,8 +700,8 @@ accept <- function(newdata, format="tibble", version = "accept3", prediction_int if (length(invalid_rows) > 0) { invalid_ids <- if ("ID" %in% colnames(newdata)) newdata$ID[invalid_rows] else invalid_rows stop(paste0("LastYrSevExacCount exceeds LastYrExacCount for patient(s) with ID(s): ", - paste(invalid_ids, collapse = ", "), - ". Severe exacerbation count cannot be greater than total exacerbation count.")) + paste(invalid_ids, collapse = ", "), + ". Severe exacerbation count cannot be greater than total exacerbation count.")) } } @@ -746,6 +746,26 @@ accept <- function(newdata, format="tibble", version = "accept3", prediction_int # Validate and normalize country code country <- toupper(country) + + # GBR has two care-setting-specific versions: + # "GBR-PRIMARY" -> primary care, ACCEPT 3.0-CPRD (accept3_cprd) + # "GBR-SPECIALTY" -> specialty care, standard accept3 GBR recalibration + # Bare "GBR" is ambiguous and defaults to specialty care with a warning. + if (country %in% c("GBR", "GBR-PRIMARY", "GBR-SPECIALTY", "GBR-SPECIALITY")) { + if (country == "GBR") { + warning("country = 'GBR' is ambiguous: defaulting to specialty care. Use 'GBR-primary' for primary care (ACCEPT 3.0-CPRD) or 'GBR-specialty' for specialty care.") + country <- "GBR-SPECIALTY" + } + + if (country == "GBR-PRIMARY") { + # Primary care: route to the CPRD-recalibrated model. + return(accept3_cprd(patientData = newdata, return_predictors = return_predictors)) + } + + # Specialty care: use the standard accept3 GBR recalibration below. + country <- "GBR" + } + supported_countries <- c("ARG", "AUS", "BRA", "CAN", "COL", "DEU", "DNK", "ESP", "FRA", "GBR", "ITA", "JPN", "KOR", "MEX", "NLD", "NOR", "SWE", "USA") if (!country %in% supported_countries && country != "XXX") { @@ -814,7 +834,7 @@ accept <- function(newdata, format="tibble", version = "accept3", prediction_int # Add predictor columns if requested if (return_predictors) { result_df <- cbind(newdata, result_df[, c("predicted_exac_probability", "predicted_exac_rate", - "predicted_severe_exac_probability", "predicted_severe_exac_rate")]) + "predicted_severe_exac_probability", "predicted_severe_exac_rate")]) } return(as_tibble(result_df)) @@ -1050,8 +1070,8 @@ predictCountProb <- function (patientResults, n = 10, shortened = TRUE){ accept3 <- function(country, ID, age, male, BMI, smoker, mMRC = NA, CVD, ICS, LABA, LAMA, LastYrExacCount, LastYrSevExacCount, FEV1, SGRQ = NA, oxygen, obs_modsev_risk) { # Create base tibble df <- tibble(country = country, ID = ID, age = age, male = male, BMI = BMI, smoker = smoker, statin = CVD, - ICS = ICS, LABA = LABA, LAMA = LAMA, LastYrExacCount = LastYrExacCount, LastYrSevExacCount = LastYrSevExacCount, - FEV1 = FEV1, oxygen = oxygen, obs_modsev_risk = obs_modsev_risk) + ICS = ICS, LABA = LABA, LAMA = LAMA, LastYrExacCount = LastYrExacCount, LastYrSevExacCount = LastYrSevExacCount, + FEV1 = FEV1, oxygen = oxygen, obs_modsev_risk = obs_modsev_risk) # Add SGRQ or mMRC depending on what's available if (!is.na(SGRQ)) { @@ -1082,3 +1102,258 @@ accept3 <- function(country, ID, age, male, BMI, smoker, mMRC = NA, CVD, ICS, LA } + + +#' Predicts COPD exacerbation risk for UK primary-care patients using ACCEPT 3.0-CPRD +#' +#' @description +#' ACCEPT 3.0-CPRD is a UK primary-care recalibration of ACCEPT 2.0 derived from +#' the CPRD (Clinical Practice Research Datalink) primary-care dataset. It applies +#' two Cox-model recalibration parameters: an optimism-corrected baseline +#' hazard (H0) and slope (beta), separately for moderate-to-severe and severe +#' exacerbations. +#' +#' When optional predictors (LABA, oxygen, ICS, LAMA, statin/CVD, BMI, smoker) +#' are missing, a CPRD-specific sequential triangular regression imputation model +#' fills them in before prediction. +#' +#' @param patientData A tibble of patients in the same format as +#' \code{samplePatients}. Mandatory columns: \code{ID}, \code{age}, +#' \code{male}, \code{FEV1}, \code{LastYrExacCount}, +#' \code{LastYrSevExacCount}, and either \code{mMRC} or \code{SGRQ}. +#' Optional columns: \code{LABA}, \code{oxygen}, \code{ICS}, \code{LAMA}, +#' \code{statin}, \code{BMI}, \code{smoker}. +#' @param return_predictors Logical. If \code{TRUE}, the input predictors are +#' returned alongside predictions. Default \code{FALSE}. +#' +#'@param quiet Logical. If \code{TRUE} (default), suppresses imputation messages. +#' Set to \code{FALSE} to report which predictors were imputed. +#' +#' @return A tibble with columns: +#' \itemize{ +#' \item \code{predicted_exac_probability} - recalibrated moderate-to-severe risk +#' \item \code{predicted_exac_rate} - corresponding rate (-log(1-p)) +#' \item \code{predicted_severe_exac_probability} - recalibrated severe risk +#' \item \code{predicted_severe_exac_rate} - corresponding rate +#' } +#' If \code{return_predictors = TRUE}, input columns are prepended. +#' +#' @details +#' ## Recalibration formula +#' For each outcome the recalibrated 1-year risk is: +#' \deqn{\hat{p}_{UK} = 1 - \exp\!\Bigl(-H_0 \cdot \exp(\beta \cdot \log(-\log(1-\hat{p}_{2})))\Bigr)} +#' where \eqn{\hat{p}_{2}} is the ACCEPT 2.0 predicted probability, and +#' \eqn{H_0} and \eqn{\beta} are the optimism-corrected parameters from the +#' CPRD Cox model (Table in manuscript): +#' \itemize{ +#' \item Moderate-to-severe: \eqn{H_0 = 0.676}, \eqn{\beta = 0.986} +#' \item Severe: \eqn{H_0 = 0.482}, \eqn{\beta = 1.124} +#' } +#' +#' ## Optional-predictor imputation +#' Missing optional predictors are imputed sequentially (triangular matrix): +#' LABA -> oxygen -> ICS -> LAMA -> statin -> BMI -> smoker. +#' Each model uses mandatory predictors (age, male, mMRC, FEV1, +#' LastYrSevExacCount, LastYrExacCount) plus any previously imputed +#' optional predictors. Binary variables are imputed with logistic regression +#' (probability rounded to 0/1); BMI with linear regression. +#' +#' @examples +#' results <- accept3_cprd(samplePatients) +#' +#' @seealso \code{\link{accept}}, \code{\link{accept2}}, \code{\link{accept3}} +#' +#' @importFrom dplyr tibble mutate select starts_with +#' @export +accept3_cprd <- function(patientData, + return_predictors = FALSE, + quiet = TRUE) { + + if (!tibble::is_tibble(patientData)) { + stop("patientData must be a tibble. Use as_tibble() to convert.") + } + if (any(patientData$LastYrSevExacCount > patientData$LastYrExacCount)) { + invalid_ids <- patientData$ID[patientData$LastYrSevExacCount > patientData$LastYrExacCount] + stop(paste0("LastYrSevExacCount exceeds LastYrExacCount for patient(s): ", + paste(invalid_ids, collapse = ", "), + ". Severe count cannot exceed total count.")) + } + + # 1. Convert mMRC -> SGRQ if needed + if (!"SGRQ" %in% colnames(patientData)) { + if ("CAT" %in% colnames(patientData)) { + warning("SGRQ not found. Using CAT score instead.") + patientData$SGRQ <- 18.87 + 1.53 * patientData$CAT + } else if ("mMRC" %in% colnames(patientData)) { + message("SGRQ not found. Using mMRC score instead.") + patientData$SGRQ <- 20.43 + 14.77 * patientData$mMRC + } else { + stop("Either mMRC, SGRQ, or CAT must be provided. None were found in patientData.") + } + } + + # 2. CPRD-specific optional-predictor imputation + # Intercepts and coefficients from Table A3 of the ACCEPT 3.0-CPRD manuscript. + # Column order in each coef vector: + # (Intercept), age, male, mMRC, FEV1, LastYrSevExacCount, + # LastYrExacCount [, LABA [, oxygen [, ICS [, LAMA [, statin [, BMI]]]]]] + + cprd_imp <- list( + + LABA = list( + binary = TRUE, + coef = c(0.418, -0.013, 0.009, 0.438, -0.004, -0.063, 0.333) + ), + + oxygen = list( + binary = TRUE, + coef = c(-7.600, 0.006, -0.197, 0.942, -0.005, 0.188, 0.040, + 0.556) # +LABA + ), + + ICS = list( + binary = TRUE, + coef = c(-1.430, -0.004, -0.157, 0.085, -0.002, -0.030, 0.168, + 3.372, 0.409) # +LABA, oxygen + ), + + LAMA = list( + binary = TRUE, + coef = c(-0.586, -0.007, 0.089, 0.337, -0.003, 0.000, 0.153, + 1.545, 0.203, -0.632) # +LABA, oxygen, ICS + ), + + statin = list( + binary = TRUE, + coef = c(-2.711, 0.015, 0.200, -0.004, 0.001, -0.032, 0.009, + 0.180, -0.196, -0.259, 0.207) # +LABA, oxygen, ICS, LAMA + ), + + BMI = list( + binary = FALSE, + clamp_low = 10, + clamp_hi = 60, + coef = c(28.944, -0.074, 0.195, 0.421, 0.027, -0.527, 0.010, + 0.451, -0.093, -0.051, -0.105, 1.629) # +LABA, oxygen, ICS, LAMA, statin + ), + + smoker = list( + binary = TRUE, + coef = c(5.449, -0.061, -0.099, 0.123, 0.000, 0.039, -0.017, + -0.056, -0.648, -0.246, 0.069, 0.031, -0.058) # all optional + ) + ) + + # Sequential imputation order matches the triangular matrix + optional_order <- c("LABA", "oxygen", "ICS", "LAMA", "statin", "BMI", "smoker") + + # Use mMRC if available, otherwise back-transform from SGRQ + if ("SGRQ" %in% colnames(patientData)) { + if (!"mMRC" %in% colnames(patientData)) { + patientData$mMRC <- (patientData$SGRQ - 20.43) / 14.77 + } else { + na_idx <- is.na(patientData$mMRC) + if (any(na_idx)) { + patientData$mMRC[na_idx] <- (patientData$SGRQ[na_idx] - 20.43) / 14.77 + } + } + } + mandatory_preds <- c("age", "male", "mMRC", "FEV1", + "LastYrSevExacCount", "LastYrExacCount") + + # Mandatory predictors drive both the optional-predictor imputation and the + # downstream accept2() call. Any remaining NA would silently propagate to NA + # predictions, so fail fast with a clear message. + missing_mandatory <- mandatory_preds[!mandatory_preds %in% colnames(patientData)] + if (length(missing_mandatory) > 0) { + stop(paste0("accept3_cprd: missing mandatory predictor column(s): ", + paste(missing_mandatory, collapse = ", "), ".")) + } + na_mandatory <- vapply(mandatory_preds, function(v) any(is.na(patientData[[v]])), logical(1)) + if (any(na_mandatory)) { + stop(paste0("accept3_cprd: mandatory predictor(s) contain missing values that could not be resolved: ", + paste(mandatory_preds[na_mandatory], collapse = ", "), + ". Provide complete values (mMRC can be supplied directly or via SGRQ).")) + } + + imputed_vars <- character(0) + + for (vname in optional_order) { + if (!vname %in% colnames(patientData) || + any(is.na(patientData[[vname]]))) { + + spec <- cprd_imp[[vname]] + pred_cols <- c(mandatory_preds, imputed_vars) + + X <- as.matrix(cbind(1, patientData[, pred_cols])) + X <- matrix(as.numeric(X), nrow = nrow(X), ncol = ncol(X)) + lp <- as.vector(X %*% spec$coef) + + if (spec$binary) { + pred_vals <- round(exp(lp) / (1 + exp(lp))) + } else { + pred_vals <- lp + if (!is.null(spec$clamp_low)) pred_vals <- pmax(pred_vals, spec$clamp_low) + if (!is.null(spec$clamp_hi)) pred_vals <- pmin(pred_vals, spec$clamp_hi) + } + + if (!vname %in% colnames(patientData)) { + patientData[[vname]] <- pred_vals + if (!quiet) message(paste0("accept3_cprd: '", vname, "' not found - imputed using CPRD model.")) + } else { + na_idx <- is.na(patientData[[vname]]) + patientData[[vname]][na_idx] <- pred_vals[na_idx] + if (any(na_idx) && !quiet) + message(paste0("accept3_cprd: ", sum(na_idx), " missing value(s) in '", + vname, "' imputed using CPRD model.")) + } + } + + # Always track all optional variables for subsequent models + # regardless of whether they needed imputation + if (vname %in% colnames(patientData)) { + imputed_vars <- c(imputed_vars, vname) + } + } + # 3. ACCEPT 2.0 predictions + accept2_preds <- accept2(patientData = patientData) + + p2_modsev <- accept2_preds$predicted_exac_probability + p2_sev <- accept2_preds$predicted_severe_exac_probability + + # 4. CPRD recalibration + # Optimism-corrected parameters from CPRD bootstrap: + # Moderate-to-severe: Baseline Hazard H0 = 0.676 (95% CI 0.671-0.681) + # Slope beta = 0.986 (95% CI 0.977-0.995) + # Severe: Baseline Hazard H0 = 0.482 (95% CI 0.472-0.495) + # Slope beta = 1.124 (95% CI 1.107-1.142) + + H0_modsev <- 0.676 + beta_modsev <- 0.986 + + H0_sev <- 0.482 + beta_sev <- 1.124 + + lp_modsev <- log(-log(1 - p2_modsev)) + lp_sev <- log(-log(1 - p2_sev)) + + cprd_modsev <- 1 - exp(-H0_modsev * exp(beta_modsev * lp_modsev)) + cprd_sev <- 1 - exp(-H0_sev * exp(beta_sev * lp_sev)) + + # 4. Assemble output + out <- dplyr::tibble( + predicted_exac_probability = round(cprd_modsev, 4), + predicted_exac_rate = round(-log(1 - cprd_modsev), 4), + predicted_severe_exac_probability = round(cprd_sev, 4), + predicted_severe_exac_rate = round(-log(1 - cprd_sev), 4) + ) + + if (return_predictors) { + out <- dplyr::bind_cols(patientData, out) + } else { + out$ID <- patientData$ID + out <- dplyr::select(out, "ID", dplyr::everything()) + } + + return(out) +} diff --git a/man/accept.Rd b/man/accept.Rd index d8db0fc..1ccc365 100644 --- a/man/accept.Rd +++ b/man/accept.Rd @@ -26,7 +26,7 @@ accept( \item{return_predictors}{default is FALSE. IF set to TRUE, returns the predictors along with prediction results.} -\item{country}{Required for accept3 version. Three-letter ISO country code (e.g., "CAN", "USA", "GBR"). Supported countries: ARG, AUS, BRA, CAN, COL, DEU, DNK, ESP, FRA, GBR, ITA, JPN, KOR, MEX, NLD, NOR, SWE, USA. For unsupported countries, provide obs_modsev_risk parameter or add it as a column in the data.} +\item{country}{Required for accept3 version. Three-letter ISO country code (e.g., "CAN", "USA", "GBR"). Supported countries: ARG, AUS, BRA, CAN, COL, DEU, DNK, ESP, FRA, GBR, ITA, JPN, KOR, MEX, NLD, NOR, SWE, USA. For the UK, two care-setting-specific versions are available: "GBR-primary" (primary care, ACCEPT 3.0-CPRD) and "GBR-specialty" (specialty care). Bare "GBR" defaults to specialty care with a warning. For unsupported countries, provide obs_modsev_risk parameter or add it as a column in the data.} \item{obs_modsev_risk}{Observed moderate-to-severe exacerbation risk for unsupported countries. Can be provided as a parameter or as a column in newdata.} diff --git a/man/accept3_cprd.Rd b/man/accept3_cprd.Rd new file mode 100644 index 0000000..4c12963 --- /dev/null +++ b/man/accept3_cprd.Rd @@ -0,0 +1,70 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/predict.R +\name{accept3_cprd} +\alias{accept3_cprd} +\title{Predicts COPD exacerbation risk for UK primary-care patients using ACCEPT 3.0-CPRD} +\usage{ +accept3_cprd(patientData, return_predictors = FALSE, quiet = TRUE) +} +\arguments{ +\item{patientData}{A tibble of patients in the same format as +\code{samplePatients}. Mandatory columns: \code{ID}, \code{age}, +\code{male}, \code{FEV1}, \code{LastYrExacCount}, +\code{LastYrSevExacCount}, and either \code{mMRC} or \code{SGRQ}. +Optional columns: \code{LABA}, \code{oxygen}, \code{ICS}, \code{LAMA}, +\code{statin}, \code{BMI}, \code{smoker}.} + +\item{return_predictors}{Logical. If \code{TRUE}, the input predictors are +returned alongside predictions. Default \code{FALSE}.} + +\item{quiet}{Logical. If \code{TRUE} (default), suppresses imputation messages. +Set to \code{FALSE} to report which predictors were imputed.} +} +\value{ +A tibble with columns: + \itemize{ + \item \code{predicted_exac_probability} - recalibrated moderate-to-severe risk + \item \code{predicted_exac_rate} - corresponding rate (-log(1-p)) + \item \code{predicted_severe_exac_probability} - recalibrated severe risk + \item \code{predicted_severe_exac_rate} - corresponding rate + } + If \code{return_predictors = TRUE}, input columns are prepended. +} +\description{ +ACCEPT 3.0-CPRD is a UK primary-care recalibration of ACCEPT 2.0 derived from +the CPRD (Clinical Practice Research Datalink) primary-care dataset. It applies +two Cox-model recalibration parameters: an optimism-corrected baseline +hazard (H0) and slope (beta), separately for moderate-to-severe and severe +exacerbations. + +When optional predictors (LABA, oxygen, ICS, LAMA, statin/CVD, BMI, smoker) +are missing, a CPRD-specific sequential triangular regression imputation model +fills them in before prediction. +} +\details{ +## Recalibration formula +For each outcome the recalibrated 1-year risk is: +\deqn{\hat{p}_{UK} = 1 - \exp\!\Bigl(-H_0 \cdot \exp(\beta \cdot \log(-\log(1-\hat{p}_{2})))\Bigr)} +where \eqn{\hat{p}_{2}} is the ACCEPT 2.0 predicted probability, and +\eqn{H_0} and \eqn{\beta} are the optimism-corrected parameters from the +CPRD Cox model (Table in manuscript): +\itemize{ + \item Moderate-to-severe: \eqn{H_0 = 0.676}, \eqn{\beta = 0.986} + \item Severe: \eqn{H_0 = 0.482}, \eqn{\beta = 1.124} +} + +## Optional-predictor imputation +Missing optional predictors are imputed sequentially (triangular matrix): +LABA -> oxygen -> ICS -> LAMA -> statin -> BMI -> smoker. +Each model uses mandatory predictors (age, male, mMRC, FEV1, +LastYrSevExacCount, LastYrExacCount) plus any previously imputed +optional predictors. Binary variables are imputed with logistic regression +(probability rounded to 0/1); BMI with linear regression. +} +\examples{ +results <- accept3_cprd(samplePatients) + +} +\seealso{ +\code{\link{accept}}, \code{\link{accept2}}, \code{\link{accept3}} +} diff --git a/tests/testthat/test-predict.R b/tests/testthat/test-predict.R index a487783..5c934bd 100644 --- a/tests/testthat/test-predict.R +++ b/tests/testthat/test-predict.R @@ -21,11 +21,11 @@ test_that("accept3 default prediction works", { expect_true("predicted_severe_exac_probability" %in% names(result)) expect_true("predicted_severe_exac_rate" %in% names(result)) expect_equal(nrow(result), 2) - + # Verify rate calculation: rate = -log(1-p) expected_rate <- -log(1 - result$predicted_exac_probability) expect_equal(result$predicted_exac_rate, expected_rate, tolerance = 1e-6) - + expected_sev_rate <- -log(1 - result$predicted_severe_exac_probability) expect_equal(result$predicted_severe_exac_rate, expected_sev_rate, tolerance = 1e-6) }) @@ -35,9 +35,9 @@ test_that("accept3 requires country parameter", { }) test_that("accept1 and accept2 warn when country parameter is provided", { - expect_warning(accept(samplePatients[1,], version="accept1", country="CAN"), + expect_warning(accept(samplePatients[1,], version="accept1", country="CAN"), "not used by accept1.*ignored") - expect_warning(accept(samplePatients[1,], version="accept2", country="USA"), + expect_warning(accept(samplePatients[1,], version="accept2", country="USA"), "not used by accept2.*ignored") }) @@ -372,3 +372,94 @@ test_that("Different SGRQ values produce different results (no conversion)", { # Different SGRQ values should produce different results expect_false(isTRUE(all.equal(result_45$predicted_exac_probability, result_50$predicted_exac_probability))) }) + + + +test_that("accept3_cprd returns correct columns and shape", { + results <- accept3_cprd(samplePatients) + expect_true(tibble::is_tibble(results)) + expect_equal(nrow(results), nrow(samplePatients)) + expect_true(all(c("ID", "predicted_exac_probability", + "predicted_exac_rate", + "predicted_severe_exac_probability", + "predicted_severe_exac_rate") %in% colnames(results))) +}) + +test_that("accept3_cprd predictions are between 0 and 1", { + results <- accept3_cprd(samplePatients) + expect_true(all(results$predicted_exac_probability > 0 & results$predicted_exac_probability < 1)) + expect_true(all(results$predicted_severe_exac_probability > 0 & results$predicted_severe_exac_probability < 1)) +}) + +test_that("accept3_cprd gives expected predictions on fixed input", { + results <- accept3_cprd(samplePatients) + expect_equal(results$predicted_exac_probability[1], 0.6986, tolerance = 0.001) + expect_equal(results$predicted_severe_exac_probability[1], 0.356, tolerance = 0.001) +}) + +test_that("accept3_cprd imputes missing optional predictors", { + patients_no_optional <- samplePatients + patients_no_optional$LABA <- NULL + patients_no_optional$oxygen <- NULL + patients_no_optional$ICS <- NULL + patients_no_optional$LAMA <- NULL + patients_no_optional$statin <- NULL + patients_no_optional$BMI <- NULL + patients_no_optional$smoker <- NULL + results <- accept3_cprd(patients_no_optional, quiet = TRUE) + expect_equal(nrow(results), nrow(samplePatients)) + expect_true(all(results$predicted_exac_probability > 0 & results$predicted_exac_probability < 1)) +}) + +test_that("accept3_cprd errors on non-tibble input", { + expect_error(accept3_cprd(as.data.frame(samplePatients))) +}) + +test_that("accept3_cprd errors when severe count exceeds total", { + patients_invalid <- samplePatients + patients_invalid$LastYrSevExacCount[1] <- 99 + expect_error(accept3_cprd(patients_invalid)) +}) + +test_that("accept3_cprd errors when a mandatory predictor is NA", { + patients_na_mand <- samplePatients + patients_na_mand$FEV1[1] <- NA + expect_error(accept3_cprd(patients_na_mand), "mandatory predictor") +}) + +test_that("accept3_cprd back-fills NA mMRC from SGRQ", { + patients_mmrc_na <- samplePatients + patients_mmrc_na$mMRC <- 2 + patients_mmrc_na$mMRC[1] <- NA + patients_mmrc_na$SGRQ <- 50 + results <- accept3_cprd(patients_mmrc_na) + expect_equal(nrow(results), nrow(samplePatients)) + expect_true(all(!is.na(results$predicted_exac_probability))) +}) + +test_that("accept(country = 'GBR-primary') routes to accept3_cprd", { + via_accept <- accept(samplePatients, country = "GBR-primary") + via_cprd <- accept3_cprd(samplePatients) + expect_equal(via_accept$predicted_exac_probability, + via_cprd$predicted_exac_probability) + expect_equal(via_accept$predicted_severe_exac_probability, + via_cprd$predicted_severe_exac_probability) +}) + +test_that("accept(country = 'GBR-specialty') uses standard accept3 recalibration", { + specialty <- accept(samplePatients, country = "GBR-specialty") + primary <- accept(samplePatients, country = "GBR-primary") + expect_true(tibble::is_tibble(specialty)) + expect_equal(nrow(specialty), nrow(samplePatients)) + # The two care settings use different recalibrations and should differ. + expect_false(isTRUE(all.equal(specialty$predicted_exac_probability, + primary$predicted_exac_probability))) +}) + +test_that("accept(country = 'GBR') warns and defaults to specialty care", { + expect_warning(gbr <- accept(samplePatients, country = "GBR"), + "ambiguous") + specialty <- accept(samplePatients, country = "GBR-specialty") + expect_equal(gbr$predicted_exac_probability, + specialty$predicted_exac_probability) +})