context("distributions.cpp unit tests")
library("flexsurv")
library("numDeriv")
library("Rcpp")
module <- Rcpp::Module('distributions', PACKAGE = "hesim")
# Exponential distribution -----------------------------------------------------
test_that("exponential", {
Exponential <- module$exponential
rate <- 2
exp <- new(Exponential, rate = rate)
# pdf
expect_equal(exp$pdf(3),
dexp(3, rate = rate))
# cdf
expect_equal(exp$cdf(3),
pexp(3, rate = rate))
# quantile
expect_equal(exp$quantile(.025),
qexp(.025, rate = rate))
# hazard
expect_equal(exp$hazard(1),
flexsurv::hexp(1, rate = rate))
# cumhazard
expect_equal(exp$cumhazard(4),
flexsurv::Hexp(4, rate = rate))
# random
set.seed(101)
r1 <- exp$random()
set.seed(101)
r2 <- rexp(1, rate = rate)
expect_equal(r1, r2)
# Truncated random
r <- replicate(10, exp$trandom(1, 5))
expect_true(all(r >= 1))
expect_true(all(r <= 5))
})
# Piecewise exponential distribution -------------------------------------------
test_that("pwexp", {
PwExp <- module$piecewise_exponential
rate <- c(.8, 1.2, 4)
time <- c(0, 5, 9)
pwexp <- new(PwExp, rate = rate, time = time)
# pdf
expect_error(pwexp$pdf(3))
# cdf
expect_error(pwexp$cdf(3))
# quantile
expect_error(exp$quantile(.025))
# hazard
expect_error(exp$hazard(1))
# cumhazard
expect_error(exp$cumhazard(4))
# random
set.seed(101)
r1 <- pwexp$random()
set.seed(101)
r2 <- rpwexp(1, rate = rate, time = time)
expect_equal(r1, r2)
# Truncated random
expect_error(
pwexp$trandom(1, 5),
paste0("hesim does not currently support sampling from a piecewise ",
"exponential distribution truncated from above.")
)
r <- replicate(10, pwexp$trandom(7, Inf))
expect_true(all(r >= 7))
r <- replicate(10, pwexp$trandom(max(time) + 1, Inf))
expect_true(all(r >= max(time) + 1))
})
# Weibull distribution (AFT) ---------------------------------------------------
test_that("weibull", {
Weibull <- module$weibull
sh <- 2; sc <- 1.2
wei <- new(Weibull, shape = sh, scale = sc)
# pdf
expect_equal(wei$pdf(3),
dweibull(3, shape = sh, scale = sc))
# cdf
expect_equal(wei$cdf(2),
pweibull(2, shape = sh, scale = sc))
# quantile
expect_equal(wei$quantile(.025),
qweibull(.025, shape = sh, scale = sc))
# hazard
expect_equal(wei$hazard(4),
flexsurv::hweibull(4, shape = sh, scale = sc))
# cumhazard
expect_equal(wei$cumhazard(4),
flexsurv::Hweibull(4, shape = sh, scale = sc))
# random
set.seed(101)
r1 <- wei$random()
set.seed(101)
r2 <- rweibull(1, shape = sh, scale = sc)
expect_equal(r1, r2)
# Truncated random
r <- replicate(10, wei$trandom(0, 1))
expect_true(all(r >= 0))
expect_true(all(r <= 1))
})
# Weibull distribution (PH) ----------------------------------------------------
test_that("weibull_ph", {
WeibullPH <- module$weibull_ph
a <- exp(0.36); m <- exp(-5)
wei <- new(WeibullPH, shape = a, scale = m)
# pdf
expect_equal(wei$pdf(4),
flexsurv::dweibullPH(4, shape = a, scale = m))
# cdf
expect_equal(wei$cdf(4),
flexsurv::pweibullPH(4, shape = a, scale = m))
# quantile
expect_equal(wei$quantile(.7),
flexsurv::qweibullPH(.7, shape = a, scale = m))
# hazard
expect_equal(wei$hazard(2),
flexsurv::hweibullPH(2, shape = a, scale = m))
# cumhazard
expect_equal(wei$cumhazard(2),
flexsurv::HweibullPH(2, shape = a, scale = m))
# random
set.seed(101)
r1 <- wei$random()
set.seed(101)
r2 <- rweibullPH(1, shape = a, scale = m)
expect_equal(r1, r2)
# Truncated random
r <- replicate(10, wei$trandom(2, 11))
expect_true(all(r >= 2))
expect_true(all(r <= 11))
})
# Weibull distribution for NMA -------------------------------------------------
test_that("weibull_nma", {
WeibullNma <- module$weibull_nma
a0 <- -2.07; a1 <- .2715
wei <- new(WeibullNma, a0 = a0, a1 = a1)
# pdf
expect_equal(wei$pdf(4),
dweibullNMA(4, a0 = a0, a1 = a1))
# cdf
expect_equal(wei$cdf(4),
pweibullNMA(4, a0 = a0, a1 = a1))
# quantile
expect_equal(wei$quantile(.7),
qweibullNMA(.7, a0 = a0, a1 = a1))
# hazard
expect_equal(wei$hazard(2),
hweibullNMA(2, a0 = a0, a1 = a1))
# cumhazard
expect_equal(wei$cumhazard(2),
HweibullNMA(2, a0 = a0, a1 = a1))
# random
set.seed(101)
r1 <- wei$random()
set.seed(101)
r2 <- rweibullNMA(1, a0 = a0, a1 = a1)
expect_equal(r1, r2)
# Truncated random
r <- replicate(10, wei$trandom(0, 11))
expect_true(all(r >= 0))
expect_true(all(r <= 11))
})
# Gamma distribution -----------------------------------------------------------
test_that("gamma", {
Gamma <- module$gamma
sh <- 2; r <- 1.4
gamma <- new(Gamma, shape = sh, rate = r)
# pdf
expect_equal(gamma$pdf(3),
dgamma(3, shape = sh, rate = r))
# cdf
expect_equal(gamma$cdf(2),
pgamma(2, shape = sh, rate = r))
# quantile
expect_equal(gamma$quantile(.025),
qgamma(.025, shape = sh, rate = r))
# hazard
expect_equal(gamma$hazard(4),
flexsurv::hgamma(4, shape = sh, rate = r))
# cumhazard
expect_equal(gamma$cumhazard(4),
flexsurv::Hgamma(4, shape = sh, rate = r))
# random
set.seed(101)
r1 <- gamma$random()
set.seed(101)
r2 <- rgamma(1, shape = sh, rate = r)
expect_equal(r1, r2)
# Truncated random
r <- replicate(10, gamma$trandom(0, 11))
expect_true(all(r >= 0))
expect_true(all(r <= 11))
})
# Lognormal distribution -------------------------------------------------------
test_that("lognormal", {
Lognormal <- module$lognormal
m <- 8; s <- 2.5
lnorm <- new(Lognormal, meanlog = m, sdlog = s)
# pdf
expect_equal(lnorm$pdf(3),
dlnorm(3, meanlog = m, sdlog = s))
# cdf
expect_equal(lnorm$cdf(3),
plnorm(3, meanlog = m, sdlog = s))
# quantile
expect_equal(lnorm$quantile(.33),
qlnorm(.33, meanlog = m, sdlog = s))
# hazard
expect_equal(lnorm$hazard(4),
flexsurv::hlnorm(4, meanlog = m, sdlog = s))
# cumhazard
expect_equal(lnorm$cumhazard(4),
flexsurv::Hlnorm(4, meanlog = m, sdlog = s))
# random
set.seed(101)
r1 <- lnorm$random()
set.seed(101)
r2 <- rlnorm(1, meanlog = m, sdlog = s)
expect_equal(r1, r2)
# Truncated random
r <- replicate(10, lnorm$trandom(10, 21))
expect_true(all(r >= 10))
expect_true(all(r <= 21))
})
# Gompertz distribution --------------------------------------------------------
test_that("gompertz", {
Gompertz <- module$gompertz
test_gompertz <- function(sh, r){
gomp <- new(Gompertz, shape = sh, rate = r)
## pdf
expect_equal(gomp$pdf(2),
flexsurv::dgompertz(2, shape = sh, rate = r))
## cdf
expect_equal(gomp$cdf(1.5),
flexsurv::pgompertz(1.5, shape = sh, rate = r))
## quantile
expect_equal(gomp$quantile(.21),
flexsurv::qgompertz(.21, shape = sh, rate = r))
## hazard
expect_equal(gomp$hazard(4),
flexsurv::hgompertz(4, shape = sh, rate = r))
## cumhazard
expect_equal(gomp$cumhazard(4),
flexsurv::Hgompertz(4, shape = sh, rate = r))
## random
set.seed(101)
r1 <- gomp$random()
set.seed(101)
r2 <- flexsurv::rgompertz(1, shape = sh, rate = r)
expect_equal(r1, r2)
## Truncated random
r <- replicate(10, gomp$trandom(2, 9))
expect_true(all(r >= 2))
expect_true(all(r <= 9))
}
test_gompertz(.05, .5) # shape > 0
test_gompertz(0, .5) # shape = 0
})
# Log-logistic distribution ----------------------------------------------------
test_that("loglogistic", {
LogLogistic <- module$loglogistic
sh <- 1; sc <- .5
llogis <- new(LogLogistic, shape = sh, scale = sc)
# pdf
expect_equal(llogis$pdf(2),
flexsurv::dllogis(2, shape = sh, scale = sc))
# cdf
expect_equal(llogis$cdf(2),
flexsurv::pllogis(2, shape = sh, scale = sc))
# quantile
expect_equal(llogis$quantile(.34),
flexsurv::qllogis(.34, shape = sh, scale = sc))
# hazard
expect_equal(llogis$hazard(6),
flexsurv::hllogis(6, shape = sh, scale = sc))
# cumhazard
expect_equal(llogis$cumhazard(6),
flexsurv::Hllogis(6, shape = sh, scale = sc))
# random
set.seed(101)
r1 <- llogis$random()
set.seed(101)
r2 <- flexsurv::rllogis(1, shape = sh, scale = sc)
expect_equal(r1, r2)
## Truncated random
r <- replicate(10, llogis$trandom(2, 9))
expect_true(all(r >= 2))
expect_true(all(r <= 9))
})
# Generalized gamma distribution -----------------------------------------------
test_that("gengamma", {
GeneralizedGamma <- module$gengamma
test_gengamma <- function(m, s, q){
gengamma <- new(GeneralizedGamma, mu = m, sigma = s, Q = q)
## pdf
expect_equal(gengamma$pdf(2),
flexsurv::dgengamma(2, mu = m, sigma = s, Q = q))
## cdf
expect_equal(gengamma$cdf(3),
flexsurv::pgengamma(3, mu = m, sigma = s, Q = q))
## quantile
if (q !=0){
expect_equal(gengamma$quantile(.3),
flexsurv::qgengamma(.3, mu = m, sigma = s, Q = q))
expect_equal(gengamma$quantile(.8),
flexsurv::qgengamma(.8, mu = m, sigma = s, Q = q))
expect_equal(gengamma$quantile(gengamma$cdf(3)), 3)
} else{ # flexsurv does not seem to be correct with Q = 0.
expect_equal(gengamma$quantile(gengamma$cdf(3)), 3)
}
## hazard
expect_equal(gengamma$hazard(1.8),
flexsurv::hgengamma(1.8, mu = m, sigma = s, Q = q))
## cumhazard
expect_equal(gengamma$cumhazard(1.2),
flexsurv::Hgengamma(1.2, mu = m, sigma = s, Q = q))
if (q == 0){ # lognormal case
## random
set.seed(101)
r1 <- gengamma$random()
set.seed(101)
r2 <- flexsurv::rgengamma(1, mu = m, sigma = s, Q = q)
expect_equal(r1, r2)
}
## Truncated random
r <- replicate(10, gengamma$trandom(2, 9))
expect_true(all(r >= 2))
expect_true(all(r <= 9))
}
test_gengamma(m = 2, s = 1.5, q = -2) # Q < 0
test_gengamma(m = 5, s = 2, q = 0) # Q = 0
test_gengamma(2, 1.1, 2) # Q > 0
})
# Spline survival distribution -------------------------------------------------
basis_cube <- function(x){
return (max(0, x^3))
}
R_linear_predict <- function(t, gamma, knots, timescale){
res <- rep(NA, length(t))
for (k in 1:length(t)){
t.scaled <- switch(timescale,
log = log(t[k]),
identity = t[k])
knot.min <- knots[1]; knot.max <- knots[length(knots)]
basis <- rep(NA, length(knots))
basis[1] <- 1; basis[2] <- t.scaled
for (j in 2:(length(knots) - 1)){
lambda.j <- (knot.max - knots[j])/(knot.max - knot.min)
basis[j + 1] <- basis_cube(t.scaled - knots[j]) -
lambda.j * basis_cube(t.scaled - knot.min) -
(1 - lambda.j) * basis_cube(t.scaled - knot.max)
}
res[k] <- basis %*% gamma
}
return (res)
}
test_that("survspline", {
SurvSpline <- module$survspline
# Scale is log hazard
test_survspline1 <- function(gamma, knots, timescale){
spline <- new(SurvSpline, gamma = gamma, knots = knots,
scale = "log_hazard", timescale = timescale,
cumhaz_method = "quad", step = -1,
random_method = "invcdf")
## hazard
expect_equal(spline$hazard(2),
exp(R_linear_predict(2, gamma, knots, timescale)))
## cumhazard
R_hazard <- function(t) {
exp(R_linear_predict(t, gamma = gamma, knots = knots,
timescale = timescale))
}
R_cumhazard <- function(t){
stats::integrate(R_hazard, 0, t)$value
}
expect_equal(spline$cumhazard(2),
R_cumhazard(2),
tolerance = .001, scale = 1)
## cdf
expect_equal(spline$cdf(.8),
1 - exp(-spline$cumhazard(.8)))
## pdf
R_cdf <- function(t){
return(1 - exp(-R_cumhazard(t)))
}
expect_equal(spline$pdf(0),
0)
expect_equal(spline$pdf(.5),
(1 - R_cdf(.5)) * R_hazard(.5))
expect_equal(spline$pdf(.5),
numDeriv::grad(R_cdf, .5))
## quantile
expect_equal(spline$quantile(spline$cdf(.55)), .55, tolerance = .001, scale = 1)
## random
expect_error(spline$random(),
NA)
## truncated random
r <- replicate(10, spline$trandom(2, 2.5))
expect_true(all(r >= 2))
expect_true(all(r <= 2.5))
}
g = c(-1.2, 1.3, .07)
k = c(.19, 1.7, 6.7)
test_survspline1(gamma = g, knots = k, timescale = "identity")
test_survspline1(gamma = g, knots = k, timescale = "log")
# scale is log cummulative hazard, log odds, or normal
test_survspline2 <- function(flexsurv_scale, timescale,
gamma = NULL, knots = NULL){
if (is.null(gamma) | is.null(knots)){
spl_fit <- flexsurvspline(Surv(recyrs, censrec) ~ 1,
data = bc, k = 1, timescale = timescale,
scale = flexsurv_scale)
gamma <- spl_fit$res.t[, "est"]
knots <- spl_fit$knots
}
hesim_scale <- switch(flexsurv_scale,
hazard = "log_cumhazard",
odds = "log_cumodds",
normal = "inv_normal"
)
spline <- new(SurvSpline, gamma = gamma, knots = knots,
scale = hesim_scale, timescale = timescale,
cumhaz_method = "quad", step = -1, random_method = "invcdf")
### pdf
expect_equal(flexsurv::dsurvspline(5, gamma = gamma, knots = knots,
scale = flexsurv_scale, timescale = timescale),
spline$pdf(5))
expect_equal(spline$pdf(0),
0)
expect_equal(spline$pdf(-5),
0)
### cdf
expect_equal(flexsurv::psurvspline(5, gamma = gamma, knots = knots,
scale = flexsurv_scale, timescale = timescale),
spline$cdf(5))
expect_equal(spline$cdf(0),
0)
expect_equal(spline$cdf(-5),
0)
### quantile
expect_equal(flexsurv::qsurvspline(.8, gamma = gamma, knots = knots,
scale = flexsurv_scale, timescale = timescale),
spline$quantile(.8),
tolerance = .001, scale = 1)
expect_equal(spline$quantile(-2),
NaN)
expect_equal(spline$quantile(0),
-Inf)
expect_equal(spline$quantile(1),
Inf)
### hazard
expect_equal(flexsurv::hsurvspline(3, gamma = gamma, knots = knots,
scale = flexsurv_scale, timescale = timescale),
spline$hazard(3))
expect_equal(spline$hazard(0),
0)
expect_equal(spline$hazard(-5),
0)
### cumhazard
expect_equal(flexsurv::Hsurvspline(3, gamma = gamma, knots = knots,
scale = flexsurv_scale, timescale = timescale),
spline$cumhazard(3))
expect_equal(spline$cumhazard(0),
0)
expect_equal(spline$cumhazard(-5),
0)
### random
set.seed(12)
r1 <- rsurvspline(1, gamma = gamma, knots = knots,
scale = flexsurv_scale, timescale = timescale)
set.seed(12)
r2 <- spline$random()
expect_equal(r1, r2, tolerance = .001, scale = 1)
## truncated random
r <- replicate(10, spline$trandom(2, 5))
expect_true(all(r >= 2))
expect_true(all(r <= 5))
}
test_survspline2(flexsurv_scale = "hazard", timescale = "log")
test_survspline2(flexsurv_scale = "hazard", timescale = "identity")
test_survspline2(flexsurv_scale = "odds", timescale = "log")
test_survspline2(flexsurv_scale = "odds", timescale = "identity")
test_survspline2(flexsurv_scale = "normal", timescale = "log")
test_survspline2(flexsurv_scale = "normal", timescale = "identity",
gamma = c(-1.2, 1.3, .07), knots = c(.19, 1.7, 6.7))
})
# Survival Fractional Polynomials ----------------------------------------------
test_that("FracPoly", {
FracPoly <- module$fracpoly
test_fracpoly <- function(gamma, powers){
R_hazard <- function(t){
return(exp(c(gamma %*% t(cbind(1, bfp(t, powers))))))
}
R_cumhazard <- function(t){
stats::integrate(R_hazard, 0, t)$value
}
R_linear_predict <- function(t){
c(gamma %*% t(cbind(1, bfp(t, powers))))
}
fp <- new(FracPoly, gamma = gamma, powers = powers,
cumhaz_method = "quad", step = -1,
random_method = "invcdf")
## hazard
expect_equal(fp$hazard(4),
exp(R_linear_predict(4)))
## cumhazard
expect_equal(fp$cumhazard(2),
R_cumhazard(2),
tolerance = .001, scale = 1)
expect_equal(fp$cumhazard(0),
0)
## cdf
expect_equal(fp$cdf(2.5),
1 - exp(-R_cumhazard(2.5)),
tolerance = .001, scale = 1)
## pdf
R_cdf <- function(t){
1 - exp(-R_cumhazard(t))
}
expect_equal(fp$pdf(.8),
numDeriv::grad(R_cdf, .8),
tolerance = .001, scale = 1)
## quantile
expect_equal(fp$quantile(fp$cdf(.45)),
.45,
tol = .001, scale = 1)
## random
### Using numeric quantile when max_x_ is infinity
set.seed(101)
r1 <- fp$random()
set.seed(101)
r2 <- fp$quantile(runif(1, 0, 1))
expect_equal(r1, r2)
### Using discrete cumulative hazard otherwise
fp$max_x_ <- 30
expect_error(fp$random(), NA) # See manual tests for correctness
## truncated random
### With max_x_ != infinity
r <- replicate(10, fp$trandom(2, 5))
expect_true(all(r >= 2 & r <= 5))
### With max_x_ == infinity
fp$max_x_ <- Inf
r <- replicate(10, fp$trandom(2, 5))
expect_true(all(r >= 2 & r <= 5))
} # end test_fracpoly
test_fracpoly(gamma = c(-1.2, -.567, 1.15),
powers = c(1, 0))
test_fracpoly(gamma = c(.2, .2),
powers = 0)
test_fracpoly(gamma = c(.2, .2, .2),
powers = c(1, 1))
test_fracpoly(gamma = c(.2, .2, .5, .5, .2, .2),
powers = c(1, 1, 2, 2, 1))
# Equivalence of fractional polynomial and gompertz
rate <- 1/5; shape <- 1
fp <- new(FracPoly, gamma = c(log(rate), shape), powers = 1,
cumhaz_method = "quad", step = -1, random_method = "invcdf")
expect_equal(fp$hazard(5),
flexsurv::hgompertz(5, shape = shape, rate = rate))
expect_equal(fp$cdf(2),
flexsurv::pgompertz(2, shape = shape, rate = rate))
# Equivalence of fractional polynomial and weibull
a0 <- -2; a1 <- 1
fp <- new(FracPoly, gamma = c(a0, a1), powers = 0,
cumhaz_method = "quad", step = -1, random_method = "invcdf")
expect_equal(fp$hazard(3),
hweibullNMA(3, a0 = a0, a1 = a1))
expect_equal(fp$cdf(4),
pweibullNMA(4, a0 = a0, a1 = a1))
})
# Point mass (fixed) distribution ----------------------------------------------
test_that("fixed", {
PointMass <- module$point_mass
est <- 3
pm <- new(PointMass, est = est)
# pdf
expect_equal(pm$pdf(est), 1)
expect_equal(pm$pdf(est - 1), 0)
# cdf
expect_equal(pm$cdf(est), 1)
expect_equal(pm$cdf(est - 1), 0)
expect_equal(pm$cdf(est + 1), 1)
# quantile
expect_error(pm$quantile(.025))
# hazard
expect_error(pm$hazard(1))
# cumhazard
expect_error(pm$cumhazard(4))
# random
expect_equal(pm$random(), est)
# Truncated random
expect_equal(pm$trandom(1, Inf), est + 1)
expect_equal(pm$trandom(1, 2), 2)
expect_equal(pm$trandom(5, Inf), est + 5)
})
# Truncated normal distribution ------------------------------------------------
test_that("rtruncnorm", {
n <- 1000
mu <- 50; sigma <- 10; lower <- 25; upper <- 60
#rtruncnorm from hesim
set.seed(10)
samp1 <- replicate(n, hesim:::C_test_rtruncnorm(mu, sigma, lower, upper))
# rtruncnorm from truncnorm package
set.seed(10)
samp3 <- truncnorm::rtruncnorm(n, lower, upper, mu, sigma)
expect_equal(samp1, samp3)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.