Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multiple variances (GSF) #395

Merged
merged 2 commits into from
Aug 12, 2024
Merged

Multiple variances (GSF) #395

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4493784
progress
gowerc Jul 31, 2024
fb67f72
progress
gowerc Jul 31, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions R/LongitudinalGSF.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,10 +81,10 @@ LongitudinalGSF <- function(
Parameter(name = "lm_gsf_mu_kg", prior = mu_kg, size = "n_arms"),
Parameter(name = "lm_gsf_mu_phi", prior = mu_phi, size = "n_arms"),

Parameter(name = "lm_gsf_omega_bsld", prior = omega_bsld, size = 1),
Parameter(name = "lm_gsf_omega_ks", prior = omega_ks, size = 1),
Parameter(name = "lm_gsf_omega_kg", prior = omega_kg, size = 1),
Parameter(name = "lm_gsf_omega_phi", prior = omega_phi, size = 1),
Parameter(name = "lm_gsf_omega_bsld", prior = omega_bsld, size = "n_studies"),
Parameter(name = "lm_gsf_omega_ks", prior = omega_ks, size = "n_arms"),
Parameter(name = "lm_gsf_omega_kg", prior = omega_kg, size = "n_arms"),
Parameter(name = "lm_gsf_omega_phi", prior = omega_phi, size = "n_arms"),

Parameter(name = "lm_gsf_sigma", prior = sigma, size = 1)
)
Expand Down
44 changes: 39 additions & 5 deletions R/SimLongitudinalGSF.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,12 @@ SimLongitudinalGSF <- function(
link_growth = 0,
link_shrinkage = 0
) {

if (length(omega_b) == 1) omega_b <- rep(omega_b, length(mu_b))
if (length(omega_s) == 1) omega_s <- rep(omega_s, length(mu_s))
if (length(omega_g) == 1) omega_g <- rep(omega_g, length(mu_g))
if (length(omega_phi) == 1) omega_phi <- rep(omega_phi, length(mu_phi))

.SimLongitudinalGSF(
times = times,
sigma = sigma,
Expand Down Expand Up @@ -112,8 +118,24 @@ setValidity(
return("The parameters `mu_s`, `mu_g` and `mu_phi` must have the same length.")
}

pairs <- list(
"omega_b" = "mu_b",
"omega_s" = "mu_s",
"omega_g" = "mu_g",
"omega_phi" = "mu_phi"
)
for (i in seq_along(pairs)) {
omega <- slot(object, names(pairs)[[i]])
mu <- slot(object, pairs[[i]])
if (!(length(omega) == length(mu))) {
return(
sprintf("`%s` must be length 1 or the same length as `%s`", omega, mu)
)
}
}

len_1_pars <- c(
"sigma", "omega_b", "omega_s", "omega_g", "omega_phi",
"sigma",
"link_dsld", "link_ttg", "link_identity", "link_growth",
"link_shrinkage"
)
Expand Down Expand Up @@ -167,13 +189,25 @@ sampleSubjects.SimLongitudinalGSF <- function(object, subjects_df) {
dplyr::distinct(.data$subject, .data$arm, .data$study) |>
dplyr::mutate(study_idx = as.numeric(.data$study)) |>
dplyr::mutate(arm_idx = as.numeric(.data$arm)) |>
dplyr::mutate(psi_b = stats::rlnorm(dplyr::n(), object@mu_b[.data$study_idx], object@omega_b)) |>
dplyr::mutate(psi_s = stats::rlnorm(dplyr::n(), object@mu_s[.data$arm_idx], object@omega_s)) |>
dplyr::mutate(psi_g = stats::rlnorm(dplyr::n(), object@mu_g[.data$arm_idx], object@omega_g)) |>
dplyr::mutate(psi_b = stats::rlnorm(
dplyr::n(),
object@mu_b[.data$study_idx],
object@omega_b[.data$study_idx]
)) |>
dplyr::mutate(psi_s = stats::rlnorm(
dplyr::n(),
object@mu_s[.data$arm_idx],
object@omega_s[.data$arm_idx]
)) |>
dplyr::mutate(psi_g = stats::rlnorm(
dplyr::n(),
object@mu_g[.data$arm_idx],
object@omega_g[.data$arm_idx]
)) |>
dplyr::mutate(psi_phi_logit = stats::rnorm(
dplyr::n(),
object@mu_phi[.data$arm_idx],
object@omega_phi
object@omega_phi[.data$arm_idx]
)) |>
dplyr::mutate(psi_phi = stats::plogis(.data$psi_phi_logit))

Expand Down
24 changes: 12 additions & 12 deletions inst/stan/lm-gsf/model.stan
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,10 @@ parameters{
vector[n_arms] lm_gsf_mu_kg;
vector[n_arms] lm_gsf_mu_phi;

real<lower={{ machine_double_eps }}> lm_gsf_omega_bsld;
real<lower={{ machine_double_eps }}> lm_gsf_omega_ks;
real<lower={{ machine_double_eps }}> lm_gsf_omega_kg;
real<lower={{ machine_double_eps }}> lm_gsf_omega_phi;
vector<lower={{ machine_double_eps }}>[n_studies] lm_gsf_omega_bsld;
vector<lower={{ machine_double_eps }}>[n_arms] lm_gsf_omega_ks;
vector<lower={{ machine_double_eps }}>[n_arms] lm_gsf_omega_kg;
vector<lower={{ machine_double_eps }}>[n_arms] lm_gsf_omega_phi;

{% if centred -%}
vector<lower={{ machine_double_eps }}>[n_subjects] lm_gsf_psi_bsld;
Expand Down Expand Up @@ -45,16 +45,16 @@ transformed parameters{

{% if not centred -%}
vector<lower={{ machine_double_eps }}>[n_subjects] lm_gsf_psi_bsld = exp(
lm_gsf_mu_bsld[subject_study_index] + (lm_gsf_eta_tilde_bsld * lm_gsf_omega_bsld)
lm_gsf_mu_bsld[subject_study_index] + (lm_gsf_eta_tilde_bsld .* lm_gsf_omega_bsld[subject_study_index])
);
vector<lower={{ machine_double_eps }}>[n_subjects] lm_gsf_psi_ks = exp(
lm_gsf_mu_ks[subject_arm_index] + (lm_gsf_eta_tilde_ks * lm_gsf_omega_ks)
lm_gsf_mu_ks[subject_arm_index] + (lm_gsf_eta_tilde_ks .* lm_gsf_omega_ks[subject_arm_index])
);
vector<lower={{ machine_double_eps }}>[n_subjects] lm_gsf_psi_kg = exp(
lm_gsf_mu_kg[subject_arm_index] + (lm_gsf_eta_tilde_kg * lm_gsf_omega_kg)
lm_gsf_mu_kg[subject_arm_index] + (lm_gsf_eta_tilde_kg .* lm_gsf_omega_kg[subject_arm_index])
);
vector[n_subjects] lm_gsf_psi_phi_logit = (
lm_gsf_mu_phi[subject_arm_index] + (lm_gsf_eta_tilde_phi * lm_gsf_omega_phi)
lm_gsf_mu_phi[subject_arm_index] + (lm_gsf_eta_tilde_phi .* lm_gsf_omega_phi[subject_arm_index])
);
{%- endif -%}
vector<
Expand Down Expand Up @@ -93,10 +93,10 @@ model {
// Source - lm-gsf/model.stan
//
{% if centred %}
lm_gsf_psi_bsld ~ lognormal(lm_gsf_mu_bsld[subject_study_index], lm_gsf_omega_bsld);
lm_gsf_psi_ks ~ lognormal(lm_gsf_mu_ks[subject_arm_index], lm_gsf_omega_ks);
lm_gsf_psi_kg ~ lognormal(lm_gsf_mu_kg[subject_arm_index], lm_gsf_omega_kg);
lm_gsf_psi_phi_logit ~ normal(lm_gsf_mu_phi[subject_arm_index], lm_gsf_omega_phi);
lm_gsf_psi_bsld ~ lognormal(lm_gsf_mu_bsld[subject_study_index], lm_gsf_omega_bsld[subject_study_index]);
lm_gsf_psi_ks ~ lognormal(lm_gsf_mu_ks[subject_arm_index], lm_gsf_omega_ks[subject_arm_index]);
lm_gsf_psi_kg ~ lognormal(lm_gsf_mu_kg[subject_arm_index], lm_gsf_omega_kg[subject_arm_index]);
lm_gsf_psi_phi_logit ~ normal(lm_gsf_mu_phi[subject_arm_index], lm_gsf_omega_phi[subject_arm_index]);
{%- endif -%}
}

72 changes: 49 additions & 23 deletions tests/testthat/test-LongitudinalGSF.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,24 @@ test_that("Can load and compile growth + shrinkage links", {
test_that("Can recover known distributional parameters from a full GSF joint model", {

skip_if_not(is_full_test())
pars <- list(
sigma = 0.01,
mu_s = log(c(0.6, 0.4)),
mu_g = log(c(0.25, 0.35)),
mu_b = log(60),
mu_phi = qlogis(c(0.4, 0.6)),
omega_b = c(0.2),
omega_s = c(0.3, 0.1),
omega_g = c(0.1, 0.3),
omega_phi = c(0.3, 0.1),
link_dsld = 0.1,
link_ttg = 0.2,
link_identity = 0,
beta_cat_B = 0.5,
beta_cat_C = -0.1,
beta_cont = 0.3,
lambda = 1 / (400 / 365)
)

set.seed(7743)
jlist <- SimJointData(
Expand All @@ -89,31 +107,31 @@ test_that("Can recover known distributional parameters from a full GSF joint mod
SimGroup(140, "Arm-B", "Study-X")
),
survival = SimSurvivalExponential(
lambda = 1 / (400 / 365),
time_max = 3,
lambda = pars$lambda,
time_max = 4,
time_step = 1 / 365,
lambda_censor = 1 / 9000,
beta_cat = c(
"A" = 0,
"B" = -0.1,
"C" = 0.5
"B" = pars$beta_cat_B,
"C" = pars$beta_cat_C
),
beta_cont = 0.3
beta_cont = pars$beta_cont
),
longitudinal = SimLongitudinalGSF(
times = c(-100, -50, 0, 1, 10, 50, 100, 150, 250, 300, 400, 500, 600) / 365,
sigma = 0.01,
mu_s = log(c(0.6, 0.4)),
mu_g = log(c(0.25, 0.35)),
mu_b = log(60),
mu_phi = qlogis(c(0.4, 0.6)),
omega_b = 0.2,
omega_s = 0.2,
omega_g = 0.2,
omega_phi = 0.2,
link_dsld = 0.1,
link_ttg = 0.2,
link_identity = 0
sigma = pars$sigma,
mu_s = pars$mu_s,
mu_g = pars$mu_g,
mu_b = pars$mu_b,
mu_phi = pars$mu_phi,
omega_b = pars$omega_b,
omega_s = pars$omega_s,
omega_g = pars$omega_g,
omega_phi = pars$omega_phi,
link_dsld = pars$link_dsld,
link_ttg = pars$link_ttg,
link_identity = pars$link_identity
),
.silent = TRUE
)
Expand Down Expand Up @@ -196,21 +214,29 @@ test_that("Can recover known distributional parameters from a full GSF joint mod

dat <- summary_post(
as.CmdStanMCMC(mp),
c("lm_gsf_mu_bsld", "lm_gsf_mu_ks", "lm_gsf_mu_kg"),
TRUE
c("lm_gsf_mu_bsld", "lm_gsf_mu_ks", "lm_gsf_mu_kg", "lm_gsf_mu_phi")
)

true_values <- c(60, 0.6, 0.4, 0.25, 0.35)
true_values <- c(pars$mu_b, pars$mu_s, pars$mu_g, pars$mu_phi)
expect_true(all(dat$q01 <= true_values))
expect_true(all(dat$q99 >= true_values))
expect_true(all(dat$ess_bulk > 100))


dat <- summary_post(
as.CmdStanMCMC(mp),
c("link_dsld", "link_ttg", "sm_exp_lambda", "lm_gsf_mu_phi")
c("lm_gsf_sigma", "lm_gsf_omega_bsld", "lm_gsf_omega_kg", "lm_gsf_omega_ks", "lm_gsf_omega_phi")
)
true_values <- c(pars$sigma, pars$omega_b, pars$omega_g, pars$omega_s, pars$omega_phi)
expect_true(all(dat$q01 <= true_values))
expect_true(all(dat$q99 >= true_values))
expect_true(all(dat$ess_bulk > 100))


true_values <- c(0.1, 0.2, 1 / (1 / (400 / 365)), qlogis(c(0.4, 0.6)))
dat <- summary_post(
as.CmdStanMCMC(mp),
c("link_dsld", "link_ttg", "sm_exp_lambda", "beta_os_cov")
)
true_values <- c(pars$link_dsld, pars$link_ttg, pars$lambda, pars$beta_cat_B, pars$beta_cat_C, pars$beta_cont)
expect_true(all(dat$q01 <= true_values))
expect_true(all(dat$q99 >= true_values))
expect_true(all(dat$ess_bulk > 100))
Expand Down
24 changes: 12 additions & 12 deletions vignettes/statistical-specification.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ Accessible via `linkIdentity()`


$$\begin{align*}
y_{ij} &\sim \mathcal{N}(SLD_{ij},\ SLD_{ij}^2 \sigma^2) \\
y_{ij} &\sim N(SLD_{ij},\ SLD_{ij}^2 \sigma^2) \\
\\
SLD_{ij} &=
\begin{cases}
Expand Down Expand Up @@ -259,17 +259,17 @@ the mean of the distribution is used.

$$
\begin{align*}
y_{ij} &\sim \mathcal{N}(SLD_{ij},\ SLD_{ij}^2 \sigma^2) \\ \\
y_{ij} &\sim N(SLD_{ij},\ SLD_{ij}^2 \sigma^2) \\ \\
SLD_{ij} &=
\begin{cases}
b_i[\phi_i e^{-s_it_{ij}} + (1-\phi_i)e^{g_i t_{ij}}] & \text{if } t_{ij}\geq 0 \\
b_i e^{g_i t_{ij}} & \text{if } t_{ij}\lt 0
\end{cases}\\
\\
b_i &\sim \text{LogNormal}(\mu_{bl(i)}, \omega_b) \\
s_i &\sim \text{LogNormal}(\mu_{sk(i)}, \omega_s) \\
g_i &\sim \text{LogNormal}(\mu_{gk(i)}, \omega_g) \\
\phi_i &\sim \text{LogitNormal}(\mu_{\phi k(i)}, \omega_\phi)
b_i &\sim \text{LogNormal}(\mu_{bl(i)}, \omega_{b l(i)}) \\
s_i &\sim \text{LogNormal}(\mu_{sk(i)}, \omega_{s k(i)}) \\
g_i &\sim \text{LogNormal}(\mu_{gk(i)}, \omega_{g k(i)}) \\
\phi_i &\sim \text{LogitNormal}(\mu_{\phi k(i)}, \omega_{\phi k(i)})
\end{align*}
$$

Expand All @@ -287,16 +287,16 @@ Where:
* $k(i)$ is the treatment arm index for subject $i$
* $l(i)$ is the study index for subject $i$
* $\mu_{\theta k(i)}$ is the population mean for parameter $\theta$ in group $k(i)$
* $\omega_{\theta}$ is the population variance for parameter $\theta$.
* $\omega_{\theta k(i)}$ is the population variance for parameter $\theta$ in group $k(i)$.


If using the non-centred parameterisation then the following alternative formulation is used:
$$
\begin{align*}
b_i &= exp(\mu_{bl(i)} + \omega_b * \eta_{b i}) \\
s_i &= exp(\mu_{sk(i)} + \omega_s * \eta_{s i}) \\
g_i &= exp(\mu_{gk(i)} + \omega_g * \eta_{g i}) \\
\phi_i &= \text{logistic}(\mu_{gk(i)} + \omega_\phi * \eta_{\phi i}) \\
b_i &= exp(\mu_{bl(i)} + \omega_{b l(i)} * \eta_{b i}) \\
s_i &= exp(\mu_{sk(i)} + \omega_{s k(i)} * \eta_{s i}) \\
g_i &= exp(\mu_{gk(i)} + \omega_{g k(i)} * \eta_{g i}) \\
\phi_i &= \text{logistic}(\mu_{gk(i)} + \omega_{\phi k(i)} * \eta_{\phi i}) \\
\\
\eta_{b i} &\sim N(0, 1)\\
\eta_{s i} &\sim N(0, 1) \\
Expand Down Expand Up @@ -386,7 +386,7 @@ non-centred parameterisation and setting the "random effects" term to be 0.

$$
\begin{align*}
y_{ij} &\sim \mathcal{N}(SLD_{ij},\ SLD_{ij}^2 \sigma^2) \\ \\
y_{ij} &\sim N(SLD_{ij},\ SLD_{ij}^2 \sigma^2) \\ \\
SLD_{ij} &=
\begin{cases} b_i e^{g_i t_{ij}} &
\text{if } t_{ij} < 0, \\
Expand Down
Loading