From ce9597ac9ee4eb9808241b2ecfcc0b69033cd5f7 Mon Sep 17 00:00:00 2001 From: Craig Gower-Page Date: Fri, 1 Mar 2024 12:13:28 +0000 Subject: [PATCH] Changed parameterisation of Loglogistic model (#275) * Changed parameterisation of log-logistic model * Added OS model formulation to statistical specification vignette --- DESCRIPTION | 1 + R/SurvivalLoglogistic.R | 14 ++-- R/simulations_os.R | 10 +-- inst/WORDLIST | 1 + inst/stan/sm-loglogistic/model.stan | 13 ++- man/SurvivalLogLogistic-class.Rd | 10 +-- man/sim_os_loglogistic.Rd | 6 +- tests/testthat/_snaps/SurvivalLoglogistic.md | 10 +-- tests/testthat/test-SurvivalLoglogistic.R | 86 +++++++++++++++++++- vignettes/statistical-specification.Rmd | 36 ++++++++ 10 files changed, 154 insertions(+), 33 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 572f14aa..80a93735 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -54,6 +54,7 @@ Suggests: testthat (>= 3.0.0), scales, lme4, + flexsurv, purrr, vdiffr, prodlim diff --git a/R/SurvivalLoglogistic.R b/R/SurvivalLoglogistic.R index f872b1a7..816f1d47 100644 --- a/R/SurvivalLoglogistic.R +++ b/R/SurvivalLoglogistic.R @@ -18,14 +18,14 @@ NULL #' @rdname SurvivalLogLogistic-class #' -#' @param lambda (`Prior`)\cr for the inverse median `lambda`. -#' @param p (`Prior`)\cr for the shape parameter `p`. -#' @param beta (`Prior`)\cr for covariates coefficients `beta`. +#' @param a (`Prior`)\cr Prior distribution for the scale parameter `a`. +#' @param b (`Prior`)\cr Prior distribution for the shape parameter `b`. +#' @param beta (`Prior`)\cr Prior distribution for covariates coefficients `beta`. #' #' @export SurvivalLogLogistic <- function( - lambda = prior_lognormal(log(0.1), 5), - p = prior_gamma(2, 5), + a = prior_lognormal(log(0.1), 5), + b = prior_gamma(2, 5), beta = prior_normal(0, 5) ) { .SurvivalLogLogistic( @@ -33,8 +33,8 @@ SurvivalLogLogistic <- function( name = "Log-Logistic", stan = StanModule("sm-loglogistic/model.stan"), parameters = ParameterList( - Parameter(name = "sm_logl_lambda", prior = lambda, size = 1), - Parameter(name = "sm_logl_p", prior = p, size = 1), + Parameter(name = "sm_loglogis_a", prior = a, size = 1), + Parameter(name = "sm_loglogis_b", prior = b, size = 1), Parameter(name = "beta_os_cov", prior = beta, size = "p_os_cov_design") ) ) diff --git a/R/simulations_os.R b/R/simulations_os.R index cd74a14e..f42e4814 100644 --- a/R/simulations_os.R +++ b/R/simulations_os.R @@ -26,15 +26,15 @@ sim_os_exponential <- function(lambda) { #' Construct a Log Hazard Function for the Log-Logistic Model #' -#' @param lambda (`number`)\cr the inverse median parameter. -#' @param p (`number`)\cr the shape parameter. +#' @param a (`number`)\cr the scale parameter. +#' @param b (`number`)\cr the shape parameter. #' #' @returns A function of `time` returning the log hazard. #' @export -sim_os_loglogistic <- function(lambda, p) { +sim_os_loglogistic <- function(a, b) { function(time) { - c1 <- log(lambda) + log(p) + (p - 1) * (log(lambda) + log(time)) - c2 <- log(1 + (lambda * time)^p) + c1 <- - log(a) + log(b) + (b - 1) * (- log(a) + log(time)) + c2 <- log(1 + (time / a)^b) return(c1 - c2) } } diff --git a/inst/WORDLIST b/inst/WORDLIST index c6450d65..74bf8995 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -42,6 +42,7 @@ ie ig ij immunotherapy +parameterise initival inits lm diff --git a/inst/stan/sm-loglogistic/model.stan b/inst/stan/sm-loglogistic/model.stan index 84403fa4..e95f9a0e 100755 --- a/inst/stan/sm-loglogistic/model.stan +++ b/inst/stan/sm-loglogistic/model.stan @@ -3,24 +3,23 @@ functions { // SurvivalLogLogistic matrix log_h0(matrix time, vector pars_os) { - real lambda = pars_os[1]; - real p = pars_os[2]; + real a = pars_os[1]; + real b = pars_os[2]; matrix[rows(time), cols(time)] result; - result =log(lambda) + log(p) + (p-1)*(log(lambda)+log(time))- log(1 + (lambda .* time) .^p) ; + result = log(b) - log(a)+ (b - 1) * (log(time) - log(a)) - log( 1 + (time./a).^b) ; return result; } } - parameters { // SurvivalLogLogistic - real sm_logl_lambda; - real sm_logl_p; + real sm_loglogis_a; + real sm_loglogis_b; } transformed parameters { // SurvivalLogLogistic - vector[2] pars_os = [sm_logl_lambda, sm_logl_p]'; + vector[2] pars_os = [sm_loglogis_a, sm_loglogis_b]'; } diff --git a/man/SurvivalLogLogistic-class.Rd b/man/SurvivalLogLogistic-class.Rd index 0ebd7886..bc36373d 100644 --- a/man/SurvivalLogLogistic-class.Rd +++ b/man/SurvivalLogLogistic-class.Rd @@ -8,17 +8,17 @@ \title{\code{SurvivalLogLogistic}} \usage{ SurvivalLogLogistic( - lambda = prior_lognormal(log(0.1), 5), - p = prior_gamma(2, 5), + a = prior_lognormal(log(0.1), 5), + b = prior_gamma(2, 5), beta = prior_normal(0, 5) ) } \arguments{ -\item{lambda}{(\code{Prior})\cr for the inverse median \code{lambda}.} +\item{a}{(\code{Prior})\cr Prior distribution for the scale parameter \code{a}.} -\item{p}{(\code{Prior})\cr for the shape parameter \code{p}.} +\item{b}{(\code{Prior})\cr Prior distribution for the shape parameter \code{b}.} -\item{beta}{(\code{Prior})\cr for covariates coefficients \code{beta}.} +\item{beta}{(\code{Prior})\cr Prior distribution for covariates coefficients \code{beta}.} } \description{ This class extends the general \code{\link{SurvivalModel}} class for using the diff --git a/man/sim_os_loglogistic.Rd b/man/sim_os_loglogistic.Rd index 486e78a8..199ce60a 100644 --- a/man/sim_os_loglogistic.Rd +++ b/man/sim_os_loglogistic.Rd @@ -4,12 +4,12 @@ \alias{sim_os_loglogistic} \title{Construct a Log Hazard Function for the Log-Logistic Model} \usage{ -sim_os_loglogistic(lambda, p) +sim_os_loglogistic(a, b) } \arguments{ -\item{lambda}{(\code{number})\cr the inverse median parameter.} +\item{a}{(\code{number})\cr the scale parameter.} -\item{p}{(\code{number})\cr the shape parameter.} +\item{b}{(\code{number})\cr the shape parameter.} } \value{ A function of \code{time} returning the log hazard. diff --git a/tests/testthat/_snaps/SurvivalLoglogistic.md b/tests/testthat/_snaps/SurvivalLoglogistic.md index 3f999a12..6a76d64c 100644 --- a/tests/testthat/_snaps/SurvivalLoglogistic.md +++ b/tests/testthat/_snaps/SurvivalLoglogistic.md @@ -6,21 +6,21 @@ Output Log-Logistic Survival Model with parameters: - sm_logl_lambda ~ lognormal(mu = -2.30259, sigma = 5) - sm_logl_p ~ gamma(alpha = 2, beta = 5) + sm_loglogis_a ~ lognormal(mu = -2.30259, sigma = 5) + sm_loglogis_b ~ gamma(alpha = 2, beta = 5) beta_os_cov ~ normal(mu = 0, sigma = 5) --- Code - x <- SurvivalLogLogistic(beta = prior_gamma(3, 4), p = prior_cauchy(0, 1)) + x <- SurvivalLogLogistic(beta = prior_gamma(3, 4), b = prior_cauchy(0, 1)) print(x) Output Log-Logistic Survival Model with parameters: - sm_logl_lambda ~ lognormal(mu = -2.30259, sigma = 5) - sm_logl_p ~ cauchy(mu = 0, sigma = 1) + sm_loglogis_a ~ lognormal(mu = -2.30259, sigma = 5) + sm_loglogis_b ~ cauchy(mu = 0, sigma = 1) beta_os_cov ~ gamma(alpha = 3, beta = 4) diff --git a/tests/testthat/test-SurvivalLoglogistic.R b/tests/testthat/test-SurvivalLoglogistic.R index 41ba8159..6f219e32 100644 --- a/tests/testthat/test-SurvivalLoglogistic.R +++ b/tests/testthat/test-SurvivalLoglogistic.R @@ -1,4 +1,88 @@ + + + + + +test_that("sim_os_loglogistic() is consistant with flexsurv", { + t <- c(1, 4, 50, 200, 600) + expect_equal( + log(flexsurv::hllogis(t, scale = 400, shape = 2)), + sim_os_loglogistic(a = 400, b = 2)(t) + ) +}) + + + +test_that("SurvivalLogLogistic can recover known values", { + + true_a <- 400 + true_b <- 2 + true_beta <- c(0.5, -0.2, 0.1) + set.seed(837) + jlist <- suppressMessages({ + simulate_joint_data( + n_arm = c(300), + times = seq(1, 2000, by = 0.5), + lambda_cen = 1 / 9000, + beta_cat = c("A" = 0, "B" = true_beta[1], "C" = true_beta[2]), + beta_cont = true_beta[3], + lm_fun = sim_lm_random_slope(link_dsld = 0, slope_mu = 0), + os_fun = sim_os_loglogistic(a = true_a, b = true_b) + ) + }) + + dat_os <- jlist$os + + jm <- JointModel( + survival = SurvivalLogLogistic( + a = prior_lognormal(log(400), 1), + b = prior_lognormal(log(2), 1) + ) + ) + + jdat <- DataJoint( + subject = DataSubject( + data = dat_os, + subject = "pt", + arm = "arm", + study = "study" + ), + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont + ) + ) + + mp <- sampleStanModel( + jm, + data = jdat, + iter_sampling = 500, + iter_warmup = 400, + chains = 1, + refresh = 0, + parallel_chains = 1 + ) + + # Variables to extract (order important) + vars <- c("sm_loglogis_a", "sm_loglogis_b", "beta_os_cov") + results_summary <- mp@results$summary(vars) + + # calculate Z-scores + par_mean <- results_summary$mean + par_sd <- results_summary$sd + par_real <- c(true_a, true_b, true_beta) + z_score <- (par_real - par_mean) / par_sd + + # Ensure Z-scores are within a reasonable margin of real values + expect_true(all(abs(z_score) <= qnorm(0.99))) + +}) + + + + + test_that("Print method for SurvivalLogLogistic works as expected", { expect_snapshot({ @@ -9,7 +93,7 @@ test_that("Print method for SurvivalLogLogistic works as expected", { expect_snapshot({ x <- SurvivalLogLogistic( beta = prior_gamma(3, 4), - p = prior_cauchy(0, 1) + b = prior_cauchy(0, 1) ) print(x) }) diff --git a/vignettes/statistical-specification.Rmd b/vignettes/statistical-specification.Rmd index 29febae4..797d1390 100644 --- a/vignettes/statistical-specification.Rmd +++ b/vignettes/statistical-specification.Rmd @@ -31,6 +31,42 @@ Please note that this document is currently a work-in-progress and does not contain complete information for this package yet. +# Survival Model Specification + +This package can only be used to fit proportional hazards models of the form: + +$$ +\log(h_i(t \mid \theta, \psi_i)) = \log(h_0(t \mid \theta)) + X_i \beta + \sum_j G_j(t \mid \psi_i) +$$ + +Where: + +- $h_0(.)$ is a parametric baseline hazard function +- $t$ is the event time +- $\theta$ is a vector of parameters that parameterise the baseline hazard distribution +- $\psi_i$ is an arbitrary vector of parameters for subject $i$ specified by the longitudinal model +- $G_j(.)$ is a link function that maps $\psi_i$ to a contribution to the log-hazard function where $j$ simply indexes the given link function +- $X_i$ is the subjects covariate design matrix +- $\beta$ is the corresponding coefficients to scale the design matrix covariates contribution to the log-hazard function + + +The following sections outline the available distributions to users which can be selected for the +baseline hazard $h_0(.)$. Please note that some of these distributions do not have +the proportional-hazards property meaning that the resulting survival model corresponding to the hazard $h_i()$ will not be of +the same parametric family as the baseline distribution with the hazard $h_0(.)$. + +## Log-Logistic Distribution + +$$ +h(x \mid a, b) = \frac +{(b/a)(x/a)^{(b-1)}} +{1 + (x/a)^b} +$$ + +Where: + +- $a > 0$ is the scale parameter +- $b > 0$ is the shape parameter # Longitudinal Model Specification