Skip to content

Commit

Permalink
Changed parameterisation of Loglogistic model (#275)
Browse files Browse the repository at this point in the history
* Changed parameterisation of log-logistic model
* Added OS model formulation to statistical specification vignette
  • Loading branch information
gowerc authored Mar 1, 2024
1 parent f5611bc commit ce9597a
Show file tree
Hide file tree
Showing 10 changed files with 154 additions and 33 deletions.
1 change: 1 addition & 0 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ Suggests:
testthat (>= 3.0.0),
scales,
lme4,
flexsurv,
purrr,
vdiffr,
prodlim
Expand Down
14 changes: 7 additions & 7 deletions R/SurvivalLoglogistic.R
Original file line number Diff line number Diff line change
Expand Up @@ -18,23 +18,23 @@ 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(
SurvivalModel(
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")
)
)
Expand Down
10 changes: 5 additions & 5 deletions R/simulations_os.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
}
}
1 change: 1 addition & 0 deletions inst/WORDLIST
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ ie
ig
ij
immunotherapy
parameterise
initival
inits
lm
Expand Down
13 changes: 6 additions & 7 deletions inst/stan/sm-loglogistic/model.stan
Original file line number Diff line number Diff line change
Expand Up @@ -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<lower={{ machine_double_eps }}> sm_logl_lambda;
real<lower={{ machine_double_eps }}> sm_logl_p;
real<lower={{ machine_double_eps }}> sm_loglogis_a;
real<lower={{ machine_double_eps }}> 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]';
}

10 changes: 5 additions & 5 deletions man/SurvivalLogLogistic-class.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

6 changes: 3 additions & 3 deletions man/sim_os_loglogistic.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

10 changes: 5 additions & 5 deletions tests/testthat/_snaps/SurvivalLoglogistic.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)

86 changes: 85 additions & 1 deletion tests/testthat/test-SurvivalLoglogistic.R
Original file line number Diff line number Diff line change
@@ -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({
Expand All @@ -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)
})
Expand Down
36 changes: 36 additions & 0 deletions vignettes/statistical-specification.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit ce9597a

Please sign in to comment.