Logit-Normal parameterisation for GSF Phi

R/LongitudinalGSF.R

@@ -31,16 +31,16 @@ NULL
 #' @param mu_bsld (`Prior`)\cr for the mean baseline value `mu_bsld`.
 #' @param mu_ks (`Prior`)\cr for the mean shrinkage rate `mu_ks`.
 #' @param mu_kg (`Prior`)\cr for the mean growth rate `mu_kg`.
+#' @param mu_phi (`Prior`)\cr for the mean proportion of cells affected by the treatment `mu_phi`.
 #'
 #' @param omega_bsld (`Prior`)\cr for the baseline value standard deviation `omega_bsld`.
 #' @param omega_ks (`Prior`)\cr for the shrinkage rate standard deviation `omega_ks`.
 #' @param omega_kg (`Prior`)\cr for the growth rate standard deviation `omega_kg`.
+#' @param omega_phi (`Prior`)\cr for the standard deviation of the proportion of cells
+#' affected by the treatment `omega_phi`.
 #'
 #' @param sigma (`Prior`)\cr for the variance of the longitudinal values `sigma`.
 #'
-#' @param a_phi (`Prior`)\cr for the alpha parameter for the fraction of cells that respond to treatment.
-#' @param b_phi (`Prior`)\cr for the beta parameter for the fraction of cells that respond to treatment.
-#'
 #' @param centred (`logical`)\cr whether to use the centred parameterization.
 #'
 #' @export
@@ -49,13 +49,12 @@ LongitudinalGSF <- function(
     mu_bsld = prior_normal(log(60), 1),
     mu_ks = prior_normal(log(0.5), 1),
     mu_kg = prior_normal(log(0.3), 1),
+    mu_phi = prior_normal(qlogis(0.5), 1),
 
     omega_bsld = prior_lognormal(log(0.2), 1),
     omega_ks = prior_lognormal(log(0.2), 1),
     omega_kg = prior_lognormal(log(0.2), 1),
-
-    a_phi = prior_lognormal(log(5), 1),
-    b_phi = prior_lognormal(log(5), 1),
+    omega_phi = prior_lognormal(log(0.2), 1),
 
     sigma = prior_lognormal(log(0.1), 1),
 
@@ -71,18 +70,12 @@ LongitudinalGSF <- function(
         Parameter(name = "lm_gsf_mu_bsld", prior = mu_bsld, size = "n_studies"),
         Parameter(name = "lm_gsf_mu_ks", prior = mu_ks, size = "n_arms"),
         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_a_phi", prior = a_phi, size = "n_arms"),
-        Parameter(name = "lm_gsf_b_phi", prior = b_phi, size = "n_arms"),
-        Parameter(
-            name = "lm_gsf_psi_phi",
-            prior = prior_init_only(prior_beta(a_phi@init, b_phi@init)),
-            size = "n_subjects"
-        ),
+        Parameter(name = "lm_gsf_omega_phi", prior = omega_phi, size = 1),
 
         Parameter(name = "lm_gsf_sigma", prior = sigma, size = 1)
     )
@@ -104,13 +97,19 @@ LongitudinalGSF <- function(
                 name = "lm_gsf_psi_kg",
                 prior = prior_init_only(prior_lognormal(mu_kg@init, omega_kg@init)),
                 size = "n_subjects"
+            ),
+            Parameter(
+                name = "lm_gsf_psi_phi_logit",
+                prior = prior_init_only(prior_normal(mu_phi@init, omega_phi@init)),
+                size = "n_subjects"
             )
         )
     } else {
         list(
             Parameter(name = "lm_gsf_eta_tilde_bsld", prior = prior_std_normal(), size = "n_subjects"),
             Parameter(name = "lm_gsf_eta_tilde_ks", prior = prior_std_normal(), size = "n_subjects"),
-            Parameter(name = "lm_gsf_eta_tilde_kg", prior = prior_std_normal(), size = "n_subjects")
+            Parameter(name = "lm_gsf_eta_tilde_kg", prior = prior_std_normal(), size = "n_subjects"),
+            Parameter(name = "lm_gsf_eta_tilde_phi", prior = prior_std_normal(), size = "n_subjects")
         )
     }
     parameters <- append(parameters, parameters_extra)

R/SimLongitudinalGSF.R

@@ -8,14 +8,15 @@ NULL
 #'
 #' @param times (`numeric`)\cr the times to generate observations at.
 #' @param sigma (`number`)\cr the variance of the longitudinal values.
-#' @param mu_s (`numeric`)\cr the mean shrinkage rates for the two treatment arms.
-#' @param mu_g (`numeric`)\cr the mean growth rates for the two treatment arms.
-#' @param mu_b (`numeric`)\cr the mean baseline values for the two treatment arms.
+#' @param mu_s (`numeric`)\cr the mean shrinkage rates.
+#' @param mu_g (`numeric`)\cr the mean growth rates.
+#' @param mu_b (`numeric`)\cr the mean baseline values.
+#' @param mu_phi (`numeric`)\cr the mean proportion of cells affected by the treatment
 #' @param omega_b (`number`)\cr the baseline value standard deviation.
 #' @param omega_s (`number`)\cr the shrinkage rate standard deviation.
 #' @param omega_g (`number`)\cr the growth rate standard deviation.
-#' @param a_phi (`number`)\cr the alpha parameter for the fraction of cells that respond to treatment.
-#' @param b_phi (`number`)\cr the beta parameter for the fraction of cells that respond to treatment.
+#' @param omega_phi (`number`)\cr for the standard deviation of the proportion of cells
+#' affected by the treatment `omega_phi`.
 #' @param link_dsld (`number`)\cr the link coefficient for the derivative contribution.
 #' @param link_ttg (`number`)\cr the link coefficient for the time-to-growth contribution.
 #' @param link_identity (`number`)\cr the link coefficient for the SLD Identity contribution.
@@ -25,11 +26,11 @@ NULL
 #' @slot mu_s (`numeric`)\cr See arguments.
 #' @slot mu_g (`numeric`)\cr See arguments.
 #' @slot mu_b (`numeric`)\cr See arguments.
+#' @slot mu_phi (`numeric`)\cr See arguments.
 #' @slot omega_b (`numeric`)\cr See arguments.
 #' @slot omega_s (`numeric`)\cr See arguments.
 #' @slot omega_g (`numeric`)\cr See arguments.
-#' @slot a_phi (`numeric`)\cr See arguments.
-#' @slot b_phi (`numeric`)\cr See arguments.
+#' @slot omega_phi (`numeric`)\cr See arguments.
 #' @slot link_dsld (`numeric`)\cr See arguments.
 #' @slot link_ttg (`numeric`)\cr See arguments.
 #' @slot link_identity (`numeric`)\cr See arguments.
@@ -45,11 +46,11 @@ NULL
         mu_s = "numeric",
         mu_g = "numeric",
         mu_b = "numeric",
-        a_phi = "numeric",
-        b_phi = "numeric",
+        mu_phi = "numeric",
         omega_b = "numeric",
         omega_s = "numeric",
         omega_g = "numeric",
+        omega_phi = "numeric",
         link_dsld = "numeric",
         link_ttg = "numeric",
         link_identity = "numeric",
@@ -65,11 +66,11 @@ SimLongitudinalGSF <- function(
     mu_s = log(c(0.6, 0.4)),
     mu_g = log(c(0.25, 0.35)),
     mu_b = log(60),
-    a_phi = c(4, 6),
-    b_phi = c(4, 6),
+    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,
     link_ttg = 0,
     link_identity = 0,
@@ -81,11 +82,11 @@ SimLongitudinalGSF <- function(
         mu_s = mu_s,
         mu_g = mu_g,
         mu_b = mu_b,
-        a_phi = a_phi,
-        b_phi = b_phi,
+        mu_phi = mu_phi,
         omega_b = omega_b,
         omega_s = omega_s,
         omega_g = omega_g,
+        omega_phi = omega_phi,
         link_dsld = link_dsld,
         link_ttg = link_ttg,
         link_identity = link_identity,
@@ -100,17 +101,15 @@ setValidity(
         par_lengths <- c(
             length(object@mu_s),
             length(object@mu_g),
-            length(object@a_phi),
-            length(object@b_phi)
+            length(object@mu_phi)
         )
         if (length(unique(par_lengths)) != 1) {
-            return("The parameters `mu_s`, `mu_g`, `a_phi`, and `b_phi` must have the same length.")
+            return("The parameters `mu_s`, `mu_g` and `mu_phi` must have the same length.")
         }
 
         len_1_pars <- c(
-            "sigma", "omega_b", "omega_s", "omega_g",
-            "link_dsld", "link_ttg", "link_identity",
-            "link_growth"
+            "sigma", "omega_b", "omega_s", "omega_g", "omega_phi",
+            "link_dsld", "link_ttg", "link_identity", "link_growth"
         )
         for (par in len_1_pars) {
             if (length(slot(object, par)) != 1) {
@@ -152,7 +151,9 @@ sampleSubjects.SimLongitudinalGSF <- function(object, subjects_df) {
         is.factor(subjects_df$study),
         is.factor(subjects_df$arm),
         length(levels(subjects_df$study)) == length(object@mu_b),
-        length(levels(subjects_df$arm)) == length(object@mu_s)
+        length(levels(subjects_df$arm)) == length(object@mu_s),
+        length(levels(subjects_df$arm)) == length(object@mu_g),
+        length(levels(subjects_df$arm)) == length(object@mu_phi)
     )
 
     res <- subjects_df |>
@@ -162,7 +163,12 @@ sampleSubjects.SimLongitudinalGSF <- function(object, subjects_df) {
         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_phi = stats::rbeta(dplyr::n(), object@a_phi[.data$arm_idx], object@b_phi[.data$arm_idx]))
+        dplyr::mutate(psi_phi_logit = stats::rnorm(
+            dplyr::n(),
+            object@mu_phi[.data$arm_idx],
+            object@omega_phi
+        )) |>
+        dplyr::mutate(psi_phi = stats::plogis(psi_phi_logit))
 
     res[, c("subject", "arm", "study", "psi_b", "psi_s", "psi_g", "psi_phi")]
 }

inst/WORDLIST

@@ -141,3 +141,4 @@ groupwise
 du
 int
 pk
+LogitNormal

inst/stan/lm-gsf/model.stan

@@ -10,26 +10,25 @@ parameters{
     vector[n_studies] lm_gsf_mu_bsld;
     vector[n_arms] lm_gsf_mu_ks;
     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;
 
 {% if centred -%}
     vector<lower={{ machine_double_eps }}>[n_subjects] lm_gsf_psi_bsld;
     vector<lower={{ machine_double_eps }}>[n_subjects] lm_gsf_psi_ks;
     vector<lower={{ machine_double_eps }}>[n_subjects] lm_gsf_psi_kg;
+    vector[n_subjects] lm_gsf_psi_phi_logit;
 {% else -%}
     vector[n_subjects] lm_gsf_eta_tilde_bsld;
     vector[n_subjects] lm_gsf_eta_tilde_ks;
     vector[n_subjects] lm_gsf_eta_tilde_kg;
+    vector[n_subjects] lm_gsf_eta_tilde_phi;
 {%- endif -%}
 
-    // Phi Parameters
-    vector<lower={{ machine_double_eps }}, upper={{ 1 - machine_double_eps }}>[n_subjects] lm_gsf_psi_phi;
-    vector<lower={{ machine_double_eps }}>[n_arms] lm_gsf_a_phi;
-    vector<lower={{ machine_double_eps }}>[n_arms] lm_gsf_b_phi;
-
     // Standard deviation of the error term
     real<lower={{ machine_double_eps }}> lm_gsf_sigma;
 
@@ -54,7 +53,14 @@ transformed parameters{
     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)
     );
+    vector[n_subjects] lm_gsf_psi_phi_logit = (
+        lm_gsf_mu_phi[subject_arm_index] + (lm_gsf_eta_tilde_phi * lm_gsf_omega_phi)
+    );
 {%- endif -%}
+    vector<
+        lower={{ machine_double_eps }},
+        upper={{ 1 - machine_double_eps }}
+    >[n_subjects] lm_gsf_psi_phi = inv_logit(lm_gsf_psi_phi_logit);
 
     vector[n_tumour_all] Ypred;
 
@@ -90,7 +96,7 @@ model {
     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);
 {%- endif -%}
-    lm_gsf_psi_phi ~ beta(lm_gsf_a_phi[subject_arm_index], lm_gsf_b_phi[subject_arm_index]);
 }
 

inst/stan/lm-gsf/quantities.stan

@@ -25,5 +25,5 @@ generated quantities {
     long_gq_pop_parameters[, 1] = exp(lm_gsf_mu_bsld[gq_long_pop_study_index]);
     long_gq_pop_parameters[, 2] = exp(lm_gsf_mu_ks[gq_long_pop_arm_index]);
     long_gq_pop_parameters[, 3] = exp(lm_gsf_mu_kg[gq_long_pop_arm_index]);
-    long_gq_pop_parameters[, 4] = lm_gsf_a_phi[gq_long_pop_arm_index] ./ (lm_gsf_a_phi[gq_long_pop_arm_index] + lm_gsf_b_phi[gq_long_pop_arm_index]);
+    long_gq_pop_parameters[, 4] = inv_logit(lm_gsf_mu_phi[gq_long_pop_arm_index]);
 }

tests/testthat/_snaps/JointModel.md

@@ -93,16 +93,16 @@
                lm_gsf_mu_bsld ~ normal(mu = 4.09434, sigma = 1)
                lm_gsf_mu_ks ~ normal(mu = -0.69315, sigma = 1)
                lm_gsf_mu_kg ~ normal(mu = -1.20397, sigma = 1)
+               lm_gsf_mu_phi ~ normal(mu = 0, sigma = 1)
                lm_gsf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1)
                lm_gsf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1)
                lm_gsf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1)
-               lm_gsf_a_phi ~ lognormal(mu = 1.60944, sigma = 1)
-               lm_gsf_b_phi ~ lognormal(mu = 1.60944, sigma = 1)
-               lm_gsf_psi_phi ~ <None>
+               lm_gsf_omega_phi ~ lognormal(mu = -1.60944, sigma = 1)
                lm_gsf_sigma ~ lognormal(mu = -2.30259, sigma = 1)
                lm_gsf_eta_tilde_bsld ~ std_normal()
                lm_gsf_eta_tilde_ks ~ std_normal()
                lm_gsf_eta_tilde_kg ~ std_normal()
+               lm_gsf_eta_tilde_phi ~ std_normal()
 
         Link:     
            No Link