Skip to content

Commit

Permalink
Clean up resmeans and test-resmeans
Browse files Browse the repository at this point in the history
  • Loading branch information
Dominic Muston committed Sep 9, 2023
1 parent 7c61591 commit e029b2c
Show file tree
Hide file tree
Showing 2 changed files with 119 additions and 89 deletions.
52 changes: 43 additions & 9 deletions R/resmeans.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,9 @@
#' # Find mean(s)
#' rmd_pf_stm(bosonc, dpam=params)
rmd_pf_stm <- function(dpam, Ty=10, starting=c(1, 0, 0)) {
# Declare local variables
Tw <- ttp.ts <- ttp.type <- ttp.spec <- NULL
ppd.ts <- ppd.type <- ppd.spec <- NULL
# Bound to aid integration in weeks
Tw <- Ty*365.25/7
# Normalize starting vector
Expand Down Expand Up @@ -109,6 +112,11 @@ prmd_pf_stm <- purrr::possibly(rmd_pf_stm, otherwise=NA_real_)
#' # Find mean(s)
#' rmd_pd_stm_cr(bosonc, dpam=params)
rmd_pd_stm_cr <- function(dpam, Ty=10, starting=c(1, 0, 0)) {
# Declare local variables
Tw <- ttp.ts <- ttp.type <- ttp.spec <- NULL
ppd.ts <- ppd.type <- ppd.spec <- NULL
pps.ts <- pps.type <- pps.spec <- NULL
S <- int_pf <- int_pd <- soj <- NULL
# Bound to aid integration in weeks
Tw <- Ty*365.25/7
# Normalize starting vector
Expand Down Expand Up @@ -175,6 +183,11 @@ prmd_pd_stm_cr <- purrr::possibly(rmd_pd_stm_cr, otherwise=NA_real_)
#' Find mean(s)
#' rmd_pd_stm_cf(bosonc, dpam=params)
rmd_pd_stm_cf <- function(dpam, Ty=10, starting=c(1, 0, 0)) {
# Declare local variables
Tw <- ttp.ts <- ttp.type <- ttp.spec <- NULL
ppd.ts <- ppd.type <- ppd.spec <- NULL
pps.ts <- pps.type <- pps.spec <- NULL
S <- int_pf <- int_pd <- soj <- NULL
# Bound to aid integration in weeks
Tw <- Ty*365.25/7
# Normalize starting vector
Expand Down Expand Up @@ -243,6 +256,9 @@ prmd_pd_stm_cf <- purrr::possibly(rmd_pd_stm_cf, otherwise=NA_real_)
#' # Find mean(s)
#' rmd_pf_psm(bosonc, dpam=params)
rmd_pf_psm <- function(dpam, Ty=10, starting=c(1, 0, 0)) {
# Declare local variables
Tw <- NULL
pfs.ts <- pfs.type <- pfs.spec <- NULL
# Bound to aid integration in weeks
Tw <- Ty*365.25/7
# Normalize starting vector
Expand Down Expand Up @@ -286,6 +302,8 @@ prmd_pf_psm <- purrr::possibly(rmd_pf_psm, otherwise=NA_real_)
#' # Find mean(s)
#' rmd_os_psm(bosonc, dpam=params)
rmd_os_psm <- function(dpam, Ty=10, starting=c(1, 0, 0)) {
# Declare local variables
Tw <- os.ts <- os.type <- os.spec <- NULL
# Bound to aid integration in weeks
Tw <- Ty*365.25/7
# Normalize starting vector
Expand Down Expand Up @@ -345,6 +363,22 @@ calc_allrmds <- function(simdat,
dpam,
cuttime = 0,
Ty = 10) {
# Declare local variables
ds <- ts.ppd <- fit.ppd <- ts.ttp <- fit.ttp <- NULL
ts.pfs <- fit.pfs <- ts.os <- fit.os <- NULL
ts.pps_cf <- fit.pps_cf <- ts.pps_cr <- fit.pps_cr <- NULL
ds <- chvalid <- pf_km <- os_km <- NULL
pfdat <- pfarea <- pfsurv <- osdat <- osarea <- ossurv <- NULL
adjTy <- pf_psm_post <- os_psm_post <- pf_stm_post <- pd_stmcf_post <- pd_stmcr_post <- NULL
pf_psm <- os_psm <- pf_stm <- os_stm_cf <- os_stm_cr <- NULL
pd_psm <- pd_stm_cf <- pd_stm_cr <- NULL
vec_pf <- vec_pd <- vec_os <- vec_model <- res <- cutadj <- NULL
# Return vector for two-piece adjustment
cutadj <- list(pf_area=pfarea, pf_surv=pfsurv, os_area=osarea, os_surv=ossurv)
# Calculate rest through differences
pd_psm_post <- os_psm_post-pf_psm_post
os_stmcf_post <- pf_stm_post+pd_stmcf_post
os_stmcr_post <- pf_stm_post+pd_stmcr_post
# For a bootstrap sample, refit all distributions
if (inclset[1]!=0) {
simdat <- simdat[inclset,]
Expand Down Expand Up @@ -419,19 +453,19 @@ calc_allrmds <- function(simdat,
pfdat <- tidyr::as_tibble(cbind(time=pf_km$time, surv=pf_km$surv)) |>
dplyr::mutate(
row = dplyr::row_number(),
itime = if_else(row==1, time, time-lag(time)),
incl = if_else(time<cuttime, 1, 0),
area = incl*surv*itime
itime = if_else(.data$row==1, .data$time, time-lag(.data$time)),
incl = if_else(.data$time<cuttime, 1, 0),
area = .data$incl*.data$surv*.data$itime
)
pfarea <- sum(pfdat$area)
pfsurv <- min(pfdat[pfdat$incl==1,]$surv)
# OS calculations
osdat <- tidyr::as_tibble(cbind(time=os_km$time, surv=os_km$surv)) |>
mutate(
row = row_number(),
itime = if_else(row==1, time, time-lag(time)),
incl = if_else(time<cuttime, 1, 0),
area = incl*surv*itime
itime = if_else(.data$row==1, .data$time, .data$time-lag(.data$time)),
incl = if_else(.data$time<cuttime, 1, 0),
area = .data$incl*.data$surv*.data$itime
)
osarea <- sum(osdat$area)
ossurv <- min(osdat[osdat$incl==1,]$surv)
Expand All @@ -450,9 +484,9 @@ calc_allrmds <- function(simdat,
pd_stmcf_post <- prmd_pd_stm_cf(dpam, Ty=adjTy, starting=starting)
pd_stmcr_post <- prmd_pd_stm_cr(dpam, Ty=adjTy, starting=starting)
# Calculate rest through differences
pd_psm_post <- os_psm_post-pf_psm_post
os_stmcf_post <- pf_stm_post+pd_stmcf_post
os_stmcr_post <- pf_stm_post+pd_stmcr_post
pd_psm_post <- os_psm_post - pf_psm_post
os_stmcf_post <- pf_stm_post + pd_stmcf_post
os_stmcr_post <- pf_stm_post + pd_stmcr_post
}
# Calculating overall means
pf_psm <- pfarea + pf_psm_post
Expand Down
156 changes: 76 additions & 80 deletions tests/testthat/test-resmeans.R
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
# Testing for resmeans.R
# ====================

library(flexsurv)
bosonc <- create_dummydata("flexbosms")
days_in_year <- 365.25
days_in_week <- 7
Expand Down Expand Up @@ -67,19 +66,19 @@ test_that("PSM results match expected durations", {
expect_equal(as.numeric(rmd_all$results[1,2]),
as.numeric(exp_psm_pf)
)
# expect_equal(as.numeric(rmd_pf_psm(params, Ty=Ty)),
# as.numeric(exp_psm_pf)
# )
# expect_equal(as.numeric(rmd_all$results[1,3]),
# as.numeric(exp_psm_pd)
# )
expect_equal(as.numeric(rmd_pf_psm(params, Ty=Ty)),
as.numeric(exp_psm_pf)
)
expect_equal(as.numeric(rmd_all$results[1,3]),
as.numeric(exp_psm_pd)
)
# rmd_pd_psm does not exist
# expect_equal(as.numeric(rmd_all$results[1,4]),
# as.numeric(exp_psm_os)
# )
# expect_equal(as.numeric(rmd_os_psm(params, Ty=Ty)),
# as.numeric(exp_psm_os)
# )
expect_equal(as.numeric(rmd_all$results[1,4]),
as.numeric(exp_psm_os)
)
expect_equal(as.numeric(rmd_os_psm(params, Ty=Ty)),
as.numeric(exp_psm_os)
)
})

# Check STM-CF calculations
Expand Down Expand Up @@ -108,24 +107,24 @@ S <- cbind(c(0,0),c(0, Tw),c(Tw, Tw))
exp_stmcf_pd <- SimplicialCubature::adaptIntegrateSimplex(int2, S)$integral
exp_stmcf_os <- exp_stmcf_pf + exp_stmcf_pd

#test_that("STM-CF results match expected durations", {
# expect_equal(as.numeric(rmd_all$results[2,2]),
# as.numeric(exp_stmcf_pf)
# )
# expect_equal(as.numeric(rmd_pf_stm(params, Ty=Ty)),
# as.numeric(exp_stmcf_pf)
# )
# expect_equal(as.numeric(rmd_all$results[2,3]),
# as.numeric(exp_stmcf_pd)
# )
# expect_equal(as.numeric(rmd_pd_stm_cf(params, Ty=Ty)),
# as.numeric(exp_stmcf_pd)
# )
# expect_equal(as.numeric(rmd_all$results[2,4]),
# as.numeric(exp_stmcf_os)
# )
test_that("STM-CF results match expected durations", {
expect_equal(as.numeric(rmd_all$results[2,2]),
as.numeric(exp_stmcf_pf)
)
expect_equal(as.numeric(rmd_pf_stm(params, Ty=Ty)),
as.numeric(exp_stmcf_pf)
)
expect_equal(as.numeric(rmd_all$results[2,3]),
as.numeric(exp_stmcf_pd)
)
expect_equal(as.numeric(rmd_pd_stm_cf(params, Ty=Ty)),
as.numeric(exp_stmcf_pd)
)
expect_equal(as.numeric(rmd_all$results[2,4]),
as.numeric(exp_stmcf_os)
)
# No rmd_os_stm_cf function
#})
})

# Check STM-CR calculations

Expand All @@ -142,24 +141,24 @@ int3 <- function(x) {
exp_stmcr_pd <- SimplicialCubature::adaptIntegrateSimplex(int3, S)$integral
exp_stmcr_os <- exp_stmcr_pf + exp_stmcr_pd

#test_that("STM-CR results match expected durations", {
# expect_equal(as.numeric(rmd_all$results[3,2]),
# as.numeric(exp_stmcr_pf)
# )
# expect_equal(as.numeric(rmd_pf_stm(params, Ty=Ty)),
# as.numeric(exp_stmcr_pf)
# )
# expect_equal(as.numeric(rmd_all$results[3,3]),
# as.numeric(exp_stmcr_pd)
# )
# expect_equal(as.numeric(rmd_pd_stm_cr(params, Ty=Ty)),
# as.numeric(exp_stmcr_pd)
# )
# expect_equal(as.numeric(rmd_all$results[3,4]),
# as.numeric(exp_stmcr_os)
# )
# No rmd_os_stm_cr function
#})
test_that("STM-CR results match expected durations", {
expect_equal(as.numeric(rmd_all$results[3,2]),
as.numeric(exp_stmcr_pf)
)
expect_equal(as.numeric(rmd_pf_stm(params, Ty=Ty)),
as.numeric(exp_stmcr_pf)
)
expect_equal(as.numeric(rmd_all$results[3,3]),
as.numeric(exp_stmcr_pd)
)
expect_equal(as.numeric(rmd_pd_stm_cr(params, Ty=Ty)),
as.numeric(exp_stmcr_pd)
)
expect_equal(as.numeric(rmd_all$results[3,4]),
as.numeric(exp_stmcr_os)
)
# No rmd_os_stm_cr function
})

# Expected results
Ty <- 10
Expand All @@ -173,39 +172,36 @@ exp_psm_os2 <- rmst_weibullPH(Tw,
shape=params$os$res[1,1],
scale=params$os$res[2,1]) # WeibullPH

#test_that("Intermediate results match expected durations", {
# expect_equal(rmd_pf_stm(params, Ty=10), exp_stmcf_pf2)
# expect_equal(rmd_pd_stm_cr(params, Ty=10), exp_stmcr_pd2)
# expect_equal(rmd_pd_stm_cf(params, Ty=10), exp_stmcf_pd2)
# expect_equal(rmd_pf_psm(params, Ty=10), exp_psm_pf2)
# expect_equal(rmd_os_psm(params, Ty=10), exp_psm_os2)
#})

#test_that("Safe functions produce the same values as originals", {
# expect_equal(rmd_pf_stm(params, Ty=15), prmd_pf_stm(params, Ty=15))
# expect_equal(rmd_pd_stm_cr(params, Ty=15), prmd_pd_stm_cr(params, Ty=15))
# expect_equal(rmd_pd_stm_cf(params, Ty=15), prmd_pd_stm_cf(params, Ty=15))
# expect_equal(rmd_pf_psm(params, Ty=15), prmd_pf_psm(params, Ty=15))
# expect_equal(rmd_os_psm(params, Ty=15), prmd_os_psm(params, Ty=15))
#})
test_that("Intermediate results match expected durations", {
expect_equal(rmd_pf_stm(params, Ty=10), exp_stmcf_pf2)
expect_equal(rmd_pd_stm_cr(params, Ty=10), exp_stmcr_pd2)
expect_equal(rmd_pd_stm_cf(params, Ty=10), exp_stmcf_pd2)
expect_equal(rmd_pf_psm(params, Ty=10), exp_psm_pf2)
expect_equal(rmd_os_psm(params, Ty=10), exp_psm_os2)
})

test_that("Safe functions produce the same values as originals", {
expect_equal(rmd_pf_stm(params, Ty=15), prmd_pf_stm(params, Ty=15))
expect_equal(rmd_pd_stm_cr(params, Ty=15), prmd_pd_stm_cr(params, Ty=15))
expect_equal(rmd_pd_stm_cf(params, Ty=15), prmd_pd_stm_cf(params, Ty=15))
expect_equal(rmd_pf_psm(params, Ty=15), prmd_pf_psm(params, Ty=15))
expect_equal(rmd_os_psm(params, Ty=15), prmd_os_psm(params, Ty=15))
})

Ty_lo <- 0.7
Ty_hi <- 1.1
#test_that("Shorter time horizon gives lower mean", {
# expect_lte(rmd_pf_stm(params, Ty=Ty_lo), rmd_pf_stm(params, Ty=Ty_hi))
# expect_lte(rmd_pd_stm_cr(params, Ty=Ty_lo), rmd_pd_stm_cr(params, Ty=Ty_hi))
# expect_lte(rmd_pd_stm_cf(params, Ty=Ty_lo), rmd_pd_stm_cf(params, Ty=Ty_hi))
# expect_lte(rmd_pf_psm(params, Ty=Ty_lo), rmd_pf_psm(params, Ty=Ty_hi))
# expect_lte(rmd_os_psm(params, Ty=Ty_lo), rmd_os_psm(params, Ty=Ty_hi))
#})

#test_that("Zero time horizon gives zero mean", {
# expect_equal(rmd_pf_stm(params, Ty=0), 0)
# expect_equal(rmd_pd_stm_cr(params, Ty=0), 0)
# expect_equal(rmd_pd_stm_cf(params, Ty=0), 0)
# expect_equal(rmd_pf_psm(params, Ty=0), 0)
# expect_equal(rmd_os_psm(params, Ty=0), 0)
#})


test_that("Shorter time horizon gives lower mean", {
expect_lte(rmd_pf_stm(params, Ty=Ty_lo), rmd_pf_stm(params, Ty=Ty_hi))
expect_lte(rmd_pd_stm_cr(params, Ty=Ty_lo), rmd_pd_stm_cr(params, Ty=Ty_hi))
expect_lte(rmd_pd_stm_cf(params, Ty=Ty_lo), rmd_pd_stm_cf(params, Ty=Ty_hi))
expect_lte(rmd_pf_psm(params, Ty=Ty_lo), rmd_pf_psm(params, Ty=Ty_hi))
expect_lte(rmd_os_psm(params, Ty=Ty_lo), rmd_os_psm(params, Ty=Ty_hi))
})

test_that("Zero time horizon gives zero mean", {
expect_equal(rmd_pf_stm(params, Ty=0), 0)
expect_equal(rmd_pd_stm_cr(params, Ty=0), 0)
expect_equal(rmd_pd_stm_cf(params, Ty=0), 0)
expect_equal(rmd_pf_psm(params, Ty=0), 0)
expect_equal(rmd_os_psm(params, Ty=0), 0)
})

0 comments on commit e029b2c

Please sign in to comment.