In this vignette we want to
compare the combination design implemented in crmPack and
explained in this vignette with the
implementation in the decider
package (Schroeter 2023). Please note
that the decider package is not available on CRAN,
therefore the vignette is only executing the code chunks if the package
is available on this system.
We are going to use the example as described in the
decider vignette here:
deciderlibrary(decider)
#> This is development version `0.0.0.9012` of the `decider` package:
#> DECIsion making in oncology Dose Escalation trials with logistic Regression.
#>
#> Attaching package: 'decider'
#> The following object is masked from 'package:crmPack':
#>
#> logitThis is the data from the historical Arm C:
historical_data <- list(
dose1 = c(0, 0, 0, 0, 0),
dose2 = c(2, 4, 8, 12, 16),
n.pat = c(3, 3, 3, 9, 12),
n.dlt = c(0, 0, 0, 1, 2),
trial = c("H1", "H1", "H1", "H1", "H1")
)The monotherapy dose grid for Arm A is:
The dose grid for compound 2 in Arm B is more sparse:
The overall dose grid for combination Arm B is therefore:
doses_of_interest <- rbind(
c(d1, rep(d1, times = length(d2))),
c(rep(0, length(d1)), rep(d2, each = length(d1)))
)The reference doses to be used in the models are:
We further need to specify the arms and types of the arms as follows:
The prior for the hypermeans is specified like this:
# Parameter Mean SD
prior_mu <- list(mu_a1 = c(logit(0.33), 2),
mu_b1 = c(0, 1), # standard normal
mu_a2 = c(logit(0.33), 2),
mu_b2 = c(0, 1), # standard normal
mu_eta = c(0, 1.121))The prior mean for \(\mu_{\alpha_{1}}\) is set to \(\text{logit}(0.33)\), which implies that we assume the reference dose has a prior median DLT rate of 33%.
Note that we use a normal prior here on the interaction parameter \(\eta\), thus allowing both positive and negative interactions. The standard deviation is set such that \(\exp(1.96 \cdot 1.121) \approx 9\), thus allowing for a 95% prior interval of \([1 / 9, 9]\) for the odds changes for a DLT at the reference dose. So \(1.121 = \log(9) / z_{0.975}\).
The prior for the between-trial heterogeneity parameters is specified like this:
# Parameter Mean SD
prior_tau <- list(tau_a1 = c(log(0.25), log(2) / 1.96),
tau_b1 = c(log(0.125), log(2) / 1.96),
tau_a2 = c(log(0.25), log(2) / 1.96),
tau_b2 = c(log(0.125), log(2) / 1.96),
tau_eta = c(log(0.125), log(2) / 1.96))These are all the log normal prior parameters for the corresponding \(\tau\) parameters. These are all “moderate” degrees of heterogeneity, according to Neuenschwander et al. (2014).
Then we look at the following scenario, where two cohorts of patients are available from Arm A:
scenario1 <- list(
dose1 = c(0.1, 0.2),
dose2 = c(0, 0),
n.pat = c(3, 3),
n.dlt = c(0, 1),
trial = c("A", "A")
)We note that the trial specification here needs to match
the name used in trials_of_interest above.
Now we can call the scenario function:
result1 <- scenario_jointBLRM(
data = scenario1,
historical.data = historical_data,
doses.of.interest = doses_of_interest,
dose.ref1 = dose_ref1,
dose.ref2 = dose_ref2,
trials.of.interest = trials_of_interest,
types.of.interest = types_of_interest,
prior.mu = prior_mu,
prior.tau = prior_tau
)
#> Warning: There were 3 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problemsWe can look at the results:
result1
#> $`trial-A`
#> mean sd q.2.5% q.50% q.97.5% P([0,0.16)) P([0.16,0.33))
#> 0.1+0 0.11481 0.10862 0.00311 0.08180 0.40465 0.74382 0.20132
#> 0.2+0 0.15349 0.12822 0.00674 0.11932 0.48034 0.61755 0.27756
#> 0.4+0 0.20650 0.15451 0.01320 0.17077 0.57876 0.47294 0.32416
#> 0.8+0 0.27469 0.18700 0.02252 0.23858 0.70410 0.33548 0.32323
#> 1.6+0 0.35422 0.22054 0.03473 0.32151 0.82456 0.22752 0.28506
#> 2.4+0 0.40297 0.23780 0.04250 0.37589 0.88227 0.18202 0.25454
#> 3.6+0 0.45126 0.25201 0.05073 0.43217 0.92488 0.14492 0.22422
#> 5+0 0.48911 0.26099 0.05768 0.47904 0.94944 0.12162 0.20159
#> 6+0 0.50938 0.26498 0.06197 0.50519 0.95994 0.11111 0.18886
#> P([0.33,1])
#> 0.1+0 0.05486
#> 0.2+0 0.10489
#> 0.4+0 0.20290
#> 0.8+0 0.34129
#> 1.6+0 0.48742
#> 2.4+0 0.56344
#> 3.6+0 0.63086
#> 5+0 0.67679
#> 6+0 0.70003
#>
#> $`trial-B`
#> mean sd q.2.5% q.50% q.97.5% P([0,0.16)) P([0.16,0.33))
#> 0.1+8 0.18544 0.12541 0.02810 0.15579 0.50934 0.51456 0.35936
#> 0.2+8 0.21994 0.14159 0.03622 0.18794 0.57506 0.41155 0.39328
#> 0.4+8 0.26671 0.16232 0.04717 0.23258 0.65782 0.30207 0.39983
#> 0.8+8 0.32727 0.18757 0.06078 0.29271 0.75632 0.20625 0.36512
#> 1.6+8 0.39966 0.21533 0.07546 0.36927 0.85382 0.13613 0.29896
#> 2.4+8 0.44499 0.23178 0.08225 0.42042 0.90118 0.11107 0.25681
#> 3.6+8 0.49007 0.24846 0.08395 0.47504 0.93919 0.09678 0.21638
#> 5+8 0.52473 0.26291 0.07941 0.52242 0.96202 0.09440 0.18670
#> 6+8 0.54263 0.27174 0.07343 0.54984 0.97211 0.09685 0.17261
#> 0.1+12 0.22352 0.12662 0.05273 0.19732 0.53946 0.36464 0.45563
#> 0.2+12 0.25643 0.14089 0.06162 0.22834 0.60001 0.27959 0.46157
#> 0.4+12 0.30106 0.15970 0.07373 0.27135 0.67736 0.19869 0.43447
#> 0.8+12 0.35894 0.18353 0.08752 0.32850 0.77139 0.13158 0.37150
#> 1.6+12 0.42817 0.21184 0.09887 0.40187 0.86582 0.09237 0.28636
#> 2.4+12 0.47125 0.23071 0.09850 0.45167 0.91328 0.08314 0.23861
#> 3.6+12 0.51329 0.25247 0.08844 0.50654 0.95043 0.08683 0.19682
#> 5+12 0.54421 0.27343 0.07095 0.55426 0.97248 0.09935 0.17011
#> 6+12 0.55927 0.28670 0.05852 0.58000 0.98165 0.11081 0.15660
#> P([0.33,1])
#> 0.1+8 0.12608
#> 0.2+8 0.19517
#> 0.4+8 0.29810
#> 0.8+8 0.42863
#> 1.6+8 0.56491
#> 2.4+8 0.63212
#> 3.6+8 0.68684
#> 5+8 0.71890
#> 6+8 0.73054
#> 0.1+12 0.17973
#> 0.2+12 0.25884
#> 0.4+12 0.36684
#> 0.8+12 0.49692
#> 1.6+12 0.62127
#> 2.4+12 0.67825
#> 3.6+12 0.71635
#> 5+12 0.73054
#> 6+12 0.73259For each trial of interest, the posterior toxicities previously designated to be of interest are shown.
Under the hood, the implementation works as follows:
post_tox_jointBLRM()
is called to sample from the posterior, which in turn usessampling_jointBLRM()
which then calls rstan::sampling() onstanmodels$jointBLRM
which is the constant Stan model sourced fromjointBLRM.stanSo we can compare this with the implementation in
crmPack which is based on JAGS.
crmPackNow we are going to define the same design and scenario in
crmPack.
We start with the monotherapy model for compound 1:
library(crmPack)
mono_model1 <- LogisticLogNormal(
mean = c(logit(0.33), 0),
cov = diag(c(2, 1)^2),
ref_dose = dose_ref1
)And for compound 2 the same:
mono_model2 <- LogisticLogNormal(
mean = c(logit(0.33), 0),
cov = diag(c(2, 1)^2),
ref_dose = dose_ref2
)Then we define the combination model:
combo_model <- TwoDrugsCombo(
list(
compound1 = mono_model1,
compound2 = mono_model2
),
gamma = 0, # prior mean for the interaction parameter
tau = 1 / (1.121^2) # prior precision for the interaction parameter
)We define the historical data which is already available:
hist_data_comp2 <- Data(
x = rep(historical_data$dose2, historical_data$n.pat),
y = c(
rep(0, sum(historical_data$n.pat) - sum(historical_data$n.dlt)),
rep(1, sum(historical_data$n.dlt))
),
doseGrid = historical_data$dose2
)
#> Used default patient IDs!
#> Used best guess cohort indices!We are going to use simple rules here (they are not relevant for the current scenario comparison):
my_stopping <- StoppingMinPatients(nPatients = 50)
my_increments <- IncrementsRelative(0, 2)
myNextBest <- NextBestNCRM(
target = c(0.16, 0.33),
overdose = c(0.33, 1),
max_overdose_prob = 0.25
)
my_cohort_size <- CohortSizeConst(size = 3)
my_increments_combo <- IncrementsComboOneDrugOnly()Then we define the design arms accordingly:
designArmA <- DesignArm(
"A",
design = Design(
data = Data(doseGrid = d1),
startingDose = d1[1],
model = mono_model1,
stopping = my_stopping,
increments = my_increments,
nextBest = myNextBest,
cohort_size = my_cohort_size
)
)
designArmB <- DesignArm(
"B",
design = DesignCombo(
data = DataCombo(doseGrid = list(compound1 = d1, compound2 = c(0, d2))),
startingDose = c(compound1 = d1[1], compound2 = 0),
model = combo_model,
stopping = my_stopping,
increments = my_increments_combo,
nextBest = myNextBest,
cohort_size = my_cohort_size
),
open_when = ArmMinDoseCondition("A", min_dose = d1[2])
)
designArmC <- HistoricalArm(
"C",
data = hist_data_comp2,
model = mono_model2
)Now we can define the hierarchical design:
design_hierarchical <- HierarchicalDesign(
designArmA,
designArmB,
designArmC,
exchangeable_parameters = list(
comp1_intercept = list(
A = "alpha0",
B = "alpha0[1]"
),
comp1_slope = list(
A = "alpha1",
B = "alpha1[1]"
),
comp2_intercept = list(
B = "alpha0[2]",
C = "alpha0"
),
comp2_slope = list(
B = "alpha1[2]",
C = "alpha1"
)
),
pool_correlations = list(
comp1 = c("comp1_intercept", "comp1_slope"),
comp2 = c("comp2_intercept", "comp2_slope")
),
pool_priors = list(
comp1_intercept = list(
mu = prior_mu$mu_a1,
tau = prior_tau$tau_a1
),
comp1_slope = list(
mu = prior_mu$mu_b1,
tau = prior_tau$tau_b1
),
comp2_intercept = list(
mu = prior_mu$mu_a2,
tau = prior_tau$tau_a2
),
comp2_slope = list(
mu = prior_mu$mu_b2,
tau = prior_tau$tau_b2
)
)
)Note that each entry in pool_correlations can correlate
exactly two scalar exchangeable parameter pools. In this example,
comp1 correlates the compound 1 intercept pool with the
compound 1 slope pool, and comp2 does the same for compound
2. Correlated blocks with three or more parameters are not currently
supported.
Then we define the scenario:
scenario_hierarchical <- HierarchicalData(
A = Data(
x = c(0.1, 0.1, 0.1, 0.2, 0.2, 0.2),
y = c(0, 0, 0, 0, 0, 1),
doseGrid = designArmA@design@data@doseGrid
),
B = designArmB@design@data,
C = designArmC@design@data
)
#> Used default patient IDs!
#> Used best guess cohort indices!And then we can use the scenario() function:
result1CrmPack <- scenario(
design_hierarchical,
data = scenario_hierarchical,
mcmcOptions = McmcOptions()
)We can look at the fit results:
result1CrmPack$fit
#> $A
#> dose middle lower upper
#> 1 0.1 0.1045962 0.0009066536 0.3706105
#> 2 0.2 0.1426731 0.0020952072 0.4437895
#> 3 0.4 0.1965217 0.0042915988 0.5571940
#> 4 0.8 0.2670607 0.0092694832 0.7016194
#> 5 1.6 0.3494819 0.0180510782 0.8497594
#> 6 2.4 0.3997139 0.0220089183 0.9155056
#> 7 3.6 0.4493329 0.0301736689 0.9542070
#> 8 5.0 0.4881666 0.0360998572 0.9716543
#> 9 6.0 0.5089692 0.0388616911 0.9790952
#>
#> $B
#> compound1 compound2 middle lower upper
#> 1 0.1 0 0.1088725 0.0005052444 0.4211303
#> 2 0.2 0 0.1454331 0.0012649891 0.4998587
#> 3 0.4 0 0.1964504 0.0028501720 0.6042296
#> 4 0.8 0 0.2639608 0.0065333146 0.7221940
#> 5 1.6 0 0.3449368 0.0140430890 0.8460175
#> 6 2.4 0 0.3953530 0.0201076564 0.9045628
#> 7 3.6 0 0.4456284 0.0256680723 0.9483034
#> 8 5.0 0 0.4851618 0.0323024732 0.9689770
#> 9 6.0 0 0.5063540 0.0356492331 0.9760133
#> 10 0.1 8 0.1704885 0.0147610510 0.4739717
#> 11 0.2 8 0.2045387 0.0177215182 0.5454017
#> 12 0.4 8 0.2521243 0.0213785787 0.6327907
#> 13 0.8 8 0.3152562 0.0271767182 0.7466187
#> 14 1.6 8 0.3913048 0.0371806177 0.8593261
#> 15 2.4 8 0.4387166 0.0431440586 0.9144964
#> 16 3.6 8 0.4856780 0.0475920935 0.9516104
#> 17 5.0 8 0.5218202 0.0485472750 0.9715697
#> 18 6.0 8 0.5405584 0.0468542930 0.9796603
#> 19 0.1 12 0.2157380 0.0477159093 0.5182767
#> 20 0.2 12 0.2480223 0.0535812842 0.5761309
#> 21 0.4 12 0.2931527 0.0588953714 0.6586895
#> 22 0.8 12 0.3530748 0.0671323909 0.7673051
#> 23 1.6 12 0.4252849 0.0730452491 0.8725445
#> 24 2.4 12 0.4700582 0.0743777639 0.9233111
#> 25 3.6 12 0.5136188 0.0687170274 0.9593125
#> 26 5.0 12 0.5458148 0.0567855866 0.9776764
#> 27 6.0 12 0.5616315 0.0477232985 0.9849400
#>
#> $C
#> dose middle lower upper
#> 1 2 0.02032274 8.849962e-07 0.1092641
#> 2 4 0.03271995 6.371554e-05 0.1339254
#> 3 8 0.06403411 4.036320e-03 0.1791379
#> 4 12 0.11211179 3.182815e-02 0.2460141
#> 5 16 0.18512166 4.890123e-02 0.3938226We can also check the probabilities to be in target and overdosing intervals:
result1CrmPack$next_best$A$probs
#> dose target overdose
#> [1,] 0.1 0.1969 0.0373
#> [2,] 0.2 0.2815 0.0863
#> [3,] 0.4 0.3397 0.1760
#> [4,] 0.8 0.3563 0.3145
#> [5,] 1.6 0.2836 0.4783
#> [6,] 2.4 0.2475 0.5616
#> [7,] 3.6 0.2131 0.6304
#> [8,] 5.0 0.1897 0.6794
#> [9,] 6.0 0.1791 0.7025
result1CrmPack$next_best$B$probs
#> compound1 compound2 target_prob overdose_prob not_eligible
#> 1 0.1 0 0.1860 0.0579 FALSE
#> 2 0.2 0 0.2522 0.1037 FALSE
#> 3 0.4 0 0.3020 0.1899 FALSE
#> 4 0.8 0 0.3203 0.3151 TRUE
#> 5 1.6 0 0.2820 0.4669 TRUE
#> 6 2.4 0 0.2466 0.5497 TRUE
#> 7 3.6 0 0.2190 0.6161 TRUE
#> 8 5.0 0 0.1945 0.6662 TRUE
#> 9 6.0 0 0.1818 0.6900 TRUE
#> 10 0.1 8 0.3329 0.1075 FALSE
#> 11 0.2 8 0.3771 0.1699 FALSE
#> 12 0.4 8 0.3965 0.2721 TRUE
#> 13 0.8 8 0.3700 0.4053 TRUE
#> 14 1.6 8 0.2955 0.5547 TRUE
#> 15 2.4 8 0.2435 0.6284 TRUE
#> 16 3.6 8 0.2016 0.6850 TRUE
#> 17 5.0 8 0.1741 0.7142 TRUE
#> 18 6.0 8 0.1619 0.7255 TRUE
#> 19 0.1 12 0.4510 0.1615 FALSE
#> 20 0.2 12 0.4681 0.2359 FALSE
#> 21 0.4 12 0.4381 0.3535 TRUE
#> 22 0.8 12 0.3673 0.4927 TRUE
#> 23 1.6 12 0.2707 0.6264 TRUE
#> 24 2.4 12 0.2248 0.6809 TRUE
#> 25 3.6 12 0.1880 0.7147 TRUE
#> 26 5.0 12 0.1585 0.7314 TRUE
#> 27 6.0 12 0.1442 0.7358 TRUEBased on this we can first compare the fit results.
Let’s look at the results for Arm A:
fitTrialADecider <- result1$`trial-A` |> as.data.frame()
fitTrialACrmPack <- result1CrmPack$fit$A
probsTrialACrmPack <- result1CrmPack$next_best$A$probs |> as.data.frame()
diffTrialA <- data.frame(
dose = fitTrialACrmPack$dose,
center = fitTrialADecider$mean - fitTrialACrmPack$middle,
lower = fitTrialADecider$`q.2.5%` - fitTrialACrmPack$lower,
upper = fitTrialADecider$`q.97.5%` - fitTrialACrmPack$upper,
target = fitTrialADecider$`P([0.16,0.33))` - probsTrialACrmPack$target,
overdose = fitTrialADecider$`P([0.33,1])` - probsTrialACrmPack$overdose
)
diffTrialA
#> dose center lower upper target overdose
#> 1 0.1 0.010213815 0.002203346 0.034039511 0.00442 0.01756
#> 2 0.2 0.010816900 0.004644793 0.036550529 -0.00394 0.01859
#> 3 0.4 0.009978256 0.008908401 0.021566017 -0.01554 0.02690
#> 4 0.8 0.007629334 0.013250517 0.002480626 -0.03307 0.02679
#> 5 1.6 0.004738106 0.016678922 -0.025199408 0.00146 0.00912
#> 6 2.4 0.003256072 0.020491082 -0.033235633 0.00704 0.00184
#> 7 3.6 0.001927142 0.020556331 -0.029327031 0.01112 0.00046
#> 8 5.0 0.000943448 0.021580143 -0.022214258 0.01189 -0.00261
#> 9 6.0 0.000410838 0.023108309 -0.019155190 0.00976 -0.00247And then the results for Arm B:
fitTrialBDecider <- result1$`trial-B` |> as.data.frame()
fitTrialBCrmPack <- result1CrmPack$fit$B |> dplyr::filter(compound2 > 0)
probsTrialBCrmPack <- result1CrmPack$next_best$B$probs |>
as.data.frame() |>
dplyr::filter(compound2 > 0)
diffTrialB <- data.frame(
dose1 = fitTrialBCrmPack$compound1,
dose2 = fitTrialBCrmPack$compound2,
center = fitTrialBDecider$mean - fitTrialBCrmPack$middle,
lower = fitTrialBDecider$`q.2.5%` - fitTrialBCrmPack$lower,
upper = fitTrialBDecider$`q.97.5%` - fitTrialBCrmPack$upper,
target = fitTrialBDecider$`P([0.16,0.33))` - probsTrialBCrmPack$target,
overdose = fitTrialBDecider$`P([0.33,1])` - probsTrialBCrmPack$overdose
)
diffTrialB
#> dose1 dose2 center lower upper target overdose
#> 1 0.1 8 0.0149514914 0.013338949 0.035368303 0.02646 0.01858
#> 2 0.2 8 0.0154012779 0.018498482 0.029658265 0.01618 0.02527
#> 3 0.4 8 0.0145856682 0.025791421 0.025029340 0.00333 0.02600
#> 4 0.8 8 0.0120138471 0.033603282 0.009701313 -0.00488 0.02333
#> 5 1.6 8 0.0083551583 0.038279382 -0.005506065 0.00346 0.01021
#> 6 2.4 8 0.0062734117 0.039105941 -0.013316360 0.01331 0.00372
#> 7 3.6 8 0.0043919516 0.036357907 -0.012420434 0.01478 0.00184
#> 8 5.0 8 0.0029097696 0.030862725 -0.009549749 0.01260 0.00470
#> 9 6.0 8 0.0020715510 0.026575707 -0.007550299 0.01071 0.00504
#> 10 0.1 12 0.0077819520 0.005014091 0.021183346 0.00463 0.01823
#> 11 0.2 12 0.0084076960 0.008038716 0.023879091 -0.00653 0.02294
#> 12 0.4 12 0.0079072693 0.014834629 0.018670486 -0.00363 0.01334
#> 13 0.8 12 0.0058652379 0.020387609 0.004084861 0.00420 0.00422
#> 14 1.6 12 0.0028850868 0.025824751 -0.006724497 0.01566 -0.00513
#> 15 2.4 12 0.0011918018 0.024122236 -0.010031106 0.01381 -0.00265
#> 16 3.6 12 -0.0003287996 0.019722973 -0.008882501 0.00882 0.00165
#> 17 5.0 12 -0.0016048451 0.014164413 -0.005196390 0.01161 -0.00086
#> 18 6.0 12 -0.0023614524 0.010796701 -0.003290025 0.01240 -0.00321So these differences look relatively small, and there does not seem to be any systematic bias in the differences.
Let’s compare the model code used in decider and
crmPack, in order to make sure that they really match and
implement the same priors and models:
deciderHere we have the following Stan model:
/*Stan model for joint BLRM
--------------------------------------------------------------------------------
Implements the joint BLRM as described in Neuenschwander et al., 2016,
"On the use of co-data in clinical trials".
A non-centered parametrization is implemented by obtaining
multivariate normals via multiplication with cholesky factors.
The cholesky decomposition is implemented by hand, as it is
available analytically in the required 2x2-case.
*/
functions{
/*counts mono observations based on input dose levels
Note: first input vector signals the component to be counted*/
int count_n_mono(vector dose_1, vector dose_2, int n_obs){
int res = 0;
for(i in 1:n_obs){
if(dose_1[i]>0 && dose_2[i]==0){
res+=1;
}
}
return res;
}
/*Computes permutation of input data, so that the first n_obs1 observations
are mono1, the subsequent n_obs2 observations are mono2, and the remaining
ones are combination therapy.
Returns matrix with two rows, first row is the permutation for sorting, and
second row contains the inverse permutation (to reverse sorted input to
normal order).*/
int[,] sort_idx(vector dose_1, vector dose_2,
int n_obs, int n_obs1, int n_obs2)
{
int res[2, n_obs] = rep_array(0, 2, n_obs);
//n_obs1/n_obs2 allow to compute offsets for sorting by counting
int cnt1 = 0;
int cnt2 = 0;
int cnt = 0;
//loop over input and save correct placement
for(i in 1:n_obs){
if(dose_1[i]>0 && dose_2[i]==0){
res[1, cnt1+1] = i;
res[2, i] = cnt1+1;
cnt1 += 1;
}else if(dose_1[i]==0 && dose_2[i]>0){
res[1, n_obs1 + 1 + cnt2] = i;
res[2, i] = n_obs1 + 1 + cnt2;
cnt2 += 1;
}else if(dose_1[i]>0 && dose_2[i]>0){
res[1, n_obs1 + n_obs2 + 1 + cnt] = i;
res[2, i] = n_obs1 + n_obs2 + 1 + cnt;
cnt += 1;
}
}
return res;
}
}
data{
//number of observations/cohorts
int<lower=0> n_obs;
//number of studies
int<lower=0> n_studies;
//number of patients for each cohort
int<lower=0> n[n_obs];
//number of DLTs for each cohort
int<lower=0> r[n_obs];
//study number for cohorts
int<lower=1> s[n_obs];
//indicates whether a MAP prior is computed
int<lower=0, upper=1> doMAP;
//indicates whether linear or saturating
//interaction term is used
int<lower=0, upper=1> saturating;
//reference doses
vector<lower=0>[2] dose_c;
//dose levels component 1 and 2 for each cohort
vector<lower=0>[n_obs] dose_1;
vector<lower=0>[n_obs] dose_2;
/*hyper priors
Notation and order of entries:
mu = (mu_alpha1, mu_beta1, mu_alpha2, mu_beta2, mu_eta)
tau = (tau_alpha1, tau_beta1, tau_alpha2, tau_beta2, tau_eta)
*/
//mean of hyper SD tau
vector[5] mean_tau;
//sd's of hyper SD tau
vector<lower=0>[5] sd_tau;
//mean of hyper mean mu
vector[5] mean_mu;
//mean of hyper sd mu
vector<lower=0>[5] sd_mu;
}
transformed data{
//internally generates a study without observations for MAP prior
int<lower=1> num_s = doMAP? n_studies+1 : n_studies;
//count number of mono observations
int<lower=0, upper=n_obs> n_obs1 = count_n_mono(dose_1, dose_2, n_obs);
int<lower=0, upper=n_obs> n_obs2 = count_n_mono(dose_2, dose_1, n_obs);
//compute sort indices (only done once per call to stan for efficiency)
int srt_idx[2, n_obs] = sort_idx(dose_1, dose_2, n_obs, n_obs1, n_obs2);
//sort by applying computed sorting permutation
int n_srt[n_obs] = n[srt_idx[1, 1:n_obs]];
int r_srt[n_obs] = r[srt_idx[1, 1:n_obs]];
int s_srt[n_obs] = s[srt_idx[1, 1:n_obs]];
//doses are also rescaled by reference dose after sorting
vector[n_obs] dose_1_srt = dose_1[srt_idx[1, 1:n_obs]]/dose_c[1];
vector[n_obs] dose_2_srt = dose_2[srt_idx[1, 1:n_obs]]/dose_c[2];
vector[n_obs] ldose_1_srt = log(dose_1_srt);
vector[n_obs] ldose_2_srt = log(dose_2_srt);
}
parameters{
//hyper SDs
real<lower=0> tau_1a;
real<lower=0> tau_1b;
real<lower=0> tau_2a;
real<lower=0> tau_2b;
real<lower=0> tau_eta;
//correlation coefficients
real<lower=-1, upper=1> rho12;
real<lower=-1, upper=1> rho34;
/*For non-centered parametrization:
Sample only raw standard normal variables. These are later transformed to
bivariate normals by multiplying with cholesky factor*/
//matrix for log(alpha_ij), log(beta_ij) and eta_j (for comp i, study j)
matrix[num_s, 5] log_ab_raw;
//for hyper means
real mu_raw[5];
}
transformed parameters{
real mu_1a;
real mu_1b;
real mu_2a;
real mu_2b;
real mu_eta;
matrix[num_s,5] log_ab;
vector<lower=0, upper=1>[n_obs] p_srt;
vector<lower=0, upper=1>[n_obs-n_obs1-n_obs2] p_2;
vector<lower=0, upper=1>[n_obs-n_obs1-n_obs2] p_1;
vector<lower=0, upper=1>[n_obs-n_obs1-n_obs2] p_0;
//transform raw hyper means to correct distribution
mu_1a = mean_mu[1] + sd_mu[1]*mu_raw[1];
mu_1b = mean_mu[2] + sd_mu[2]*mu_raw[2];
mu_2a = mean_mu[3] + sd_mu[3]*mu_raw[3];
mu_2b = mean_mu[4] + sd_mu[4]*mu_raw[4];
mu_eta = mean_mu[5] + sd_mu[5]*mu_raw[5];
/*Hard-coded matrix multiplication with lower cholesky factor
of covariance matrix. This can be done without saving the
cholesky factor itself, as it is available analytically.
The following means:
log_ab = mu + L*log_ab_raw,
where L is a lower triangular matrix with L*L^T=Sigma,
for a covariance matrix Sigma.
Note: For general
Sigma = tau_1^2 rho*tau_1*tau_2
rho*tau_1*tau_2 tau_2^2
the lower cholesky factor is
L = tau_1 0
tau_2*rho tau_2*squareroot(1-rho^2)
*/
log_ab[1:num_s,1] = mu_1a + tau_1a*log_ab_raw[1:num_s, 1];
log_ab[1:num_s,2] = mu_1b + tau_1b*rho12*log_ab_raw[1:num_s, 1] +
tau_1b*sqrt(1-square(rho12))*log_ab_raw[1:num_s, 2];
log_ab[1:num_s,3] = mu_2a + tau_2a*log_ab_raw[1:num_s, 3];
log_ab[1:num_s,4] = mu_2b + tau_2b*rho34*log_ab_raw[1:num_s, 3] +
tau_2b*sqrt(1-square(rho34))*log_ab_raw[1:num_s, 4];
log_ab[1:num_s,5] = mu_eta + tau_eta*log_ab_raw[1:num_s, 5];
//toxicity models for mono and combination treatment are vectorized
if(n_obs1>0){
//treatments mono 1
p_srt[1:n_obs1] = inv_logit(log_ab[s_srt[1:n_obs1],1] +
(exp(log_ab[s_srt[1:n_obs1],2]).*
ldose_1_srt[1:n_obs1]));
}
if(n_obs2>0){
//treatments mono 2
p_srt[(n_obs1+1):(n_obs1+n_obs2)] =
inv_logit(log_ab[s_srt[(n_obs1+1):(n_obs1 + n_obs2)],3] +
(exp(log_ab[s_srt[(n_obs1+1): (n_obs1 + n_obs2)],4]).*
ldose_2_srt[(n_obs1+1): (n_obs1 + n_obs2)]));
}
if(n_obs-n_obs1-n_obs2>0){
//treatments combination
p_2[1 : (n_obs-n_obs1-n_obs2)] =
inv_logit(log_ab[s_srt[(n_obs1 + n_obs2 + 1) : n_obs],3] +
(exp(log_ab[s_srt[(n_obs1 + n_obs2 + 1) : n_obs],4]).*
ldose_2_srt[(n_obs1 + n_obs2 + 1) : n_obs]));
p_1[1 : (n_obs-n_obs1-n_obs2)] =
inv_logit(log_ab[s_srt[(n_obs1 + n_obs2 + 1) : n_obs],1] +
(exp(log_ab[s_srt[(n_obs1 + n_obs2 + 1) : n_obs],2]).*
ldose_1_srt[(n_obs1 + n_obs2 + 1) : n_obs]));
p_0[1 :(n_obs-n_obs1-n_obs2)] = p_1[1 : (n_obs-n_obs1-n_obs2)] +
p_2[1 : (n_obs-n_obs1-n_obs2)] -
(p_1[1 : (n_obs-n_obs1-n_obs2)] .*
p_2[1 : (n_obs-n_obs1-n_obs2)]);
if(saturating){
p_srt[(n_obs1 + n_obs2 + 1) : n_obs] =
inv_logit(logit(p_0[1 : (n_obs-n_obs1-n_obs2)]) +
(2*log_ab[s_srt[(n_obs1 + n_obs2 + 1) : n_obs],5].*
(dose_1_srt[(n_obs1 + n_obs2 + 1) : n_obs].*
dose_2_srt[(n_obs1 + n_obs2 + 1) : n_obs] )./
(1 + dose_1_srt[(n_obs1 + n_obs2 + 1) : n_obs].*
dose_2_srt[(n_obs1 + n_obs2 + 1) : n_obs])
));
}else{
p_srt[(n_obs1 + n_obs2 + 1) : n_obs] =
inv_logit(logit(p_0[1 : (n_obs-n_obs1-n_obs2)]) +
log_ab[s_srt[(n_obs1 + n_obs2 + 1) : n_obs],5].*
dose_1_srt[(n_obs1 + n_obs2 + 1) : n_obs].*
dose_2_srt[(n_obs1 + n_obs2 + 1) : n_obs] );
}
}
}
model{
//priors for hyper means (non-centered)
mu_raw ~ std_normal();
//priors for hyper SD
tau_1a ~ lognormal(mean_tau[1], sd_tau[1]);
tau_1b ~ lognormal(mean_tau[2], sd_tau[2]);
tau_2a ~ lognormal(mean_tau[3], sd_tau[3]);
tau_2b ~ lognormal(mean_tau[4], sd_tau[4]);
tau_eta ~ lognormal(mean_tau[5], sd_tau[5]);
//priors for correlation coefficients
rho12 ~ uniform(-1,1);
rho34 ~ uniform(-1,1);
//priors for regression parameters (non-centered)
for(k in 1:num_s){
log_ab_raw[k, 1:5] ~ std_normal();
}
//binomial likelihood
r_srt ~ binomial(n_srt, p_srt);
}
generated quantities{
//just to provide the sorted toxicity parameters as output
vector<lower=0, upper=1>[n_obs] p = p_srt[srt_idx[2,1:n_obs]];
}crmPackHere we have the following JAGS model:
{
for (i in 1:nObs_A) {
logit(p_A[i]) <- alpha0_A + alpha1_A * log(x_A[i]/ref_dose_A)
y_A[i] ~ dbern(p_A[i])
}
for (i in 1:nObs_B) {
x_drug1_B[i] <- x_B[i, 1L]
}
for (i in 1:nObs_B) {
logit(p_drug1_B[i]) <- alpha0_drug1_B + alpha1_drug1_B *
log(x_drug1_B[i]/ref_dose_drug1_B)
p_single_B[i, 1L] <- p_drug1_B[i]
}
for (i in 1:nObs_B) {
x_drug2_B[i] <- x_B[i, 2L]
}
for (i in 1:nObs_B) {
logit(p_drug2_B[i]) <- alpha0_drug2_B + alpha1_drug2_B *
log(x_drug2_B[i]/ref_dose_drug2_B)
p_single_B[i, 2L] <- p_drug2_B[i]
}
for (i in 1:nObs_B) {
combo_interaction_B[i] <- x_drug1_B[i]/ref_dose_drug1_B *
(x_drug2_B[i]/ref_dose_drug2_B)
}
for (i in 1:nObs_B) {
p0_B[i] <- p_single_B[i, 1] + p_single_B[i, 2] - p_single_B[i,
1] * p_single_B[i, 2]
logit(p_B[i]) <- log(p0_B[i]/(1 - p0_B[i])) + eta_B *
combo_interaction_B[i]
y_B[i] ~ dbern(p_B[i])
}
for (i in 1:nObs_C) {
logit(p_C[i]) <- alpha0_C + alpha1_C * log(x_C[i]/ref_dose_C)
y_C[i] ~ dbern(p_C[i])
}
}
{
alpha0_A <- theta_A[1]
alpha1_A <- exp(theta_A[2])
alpha0_drug1_B <- theta_drug1_B[1]
alpha1_drug1_B <- exp(theta_drug1_B[2])
alpha0_drug2_B <- theta_drug2_B[1]
alpha1_drug2_B <- exp(theta_drug2_B[2])
alpha0_B[1L] <- alpha0_drug1_B
alpha0_B[2L] <- alpha0_drug2_B
alpha1_B[1L] <- alpha1_drug1_B
alpha1_B[2L] <- alpha1_drug2_B
eta_B ~ dnorm(eta_gamma_B, eta_tau_B)
alpha0_C <- theta_C[1]
alpha1_C <- exp(theta_C[2])
theta_A[1:2] ~ dmnorm(mu_comp1_corr[], prec_comp1_corr[,
])
theta_drug1_B[1:2] ~ dmnorm(mu_comp1_corr[], prec_comp1_corr[,
])
mu_comp1_corr[1] <- mu_comp1_intercept
mu_comp1_corr[2] <- mu_comp1_slope
rho_comp1 ~ dunif(rho_comp1_lower, rho_comp1_upper)
prec_comp1_corr[1, 1] <- 1/(pow(tau_comp1_intercept, 2) *
(1 - pow(rho_comp1, 2)))
prec_comp1_corr[2, 2] <- 1/(pow(tau_comp1_slope, 2) * (1 -
pow(rho_comp1, 2)))
prec_comp1_corr[1, 2] <- -rho_comp1/(tau_comp1_intercept *
tau_comp1_slope * (1 - pow(rho_comp1, 2)))
prec_comp1_corr[2, 1] <- prec_comp1_corr[1, 2]
mu_comp1_intercept ~ dnorm(mu_comp1_intercept_mean, pow(mu_comp1_intercept_sd,
-2))
tau_comp1_intercept ~ dlnorm(tau_comp1_intercept_meanlog,
pow(tau_comp1_intercept_sdlog, -2))
mu_comp1_slope ~ dnorm(mu_comp1_slope_mean, pow(mu_comp1_slope_sd,
-2))
tau_comp1_slope ~ dlnorm(tau_comp1_slope_meanlog, pow(tau_comp1_slope_sdlog,
-2))
theta_drug2_B[1:2] ~ dmnorm(mu_comp2_corr[], prec_comp2_corr[,
])
theta_C[1:2] ~ dmnorm(mu_comp2_corr[], prec_comp2_corr[,
])
mu_comp2_corr[1] <- mu_comp2_intercept
mu_comp2_corr[2] <- mu_comp2_slope
rho_comp2 ~ dunif(rho_comp2_lower, rho_comp2_upper)
prec_comp2_corr[1, 1] <- 1/(pow(tau_comp2_intercept, 2) *
(1 - pow(rho_comp2, 2)))
prec_comp2_corr[2, 2] <- 1/(pow(tau_comp2_slope, 2) * (1 -
pow(rho_comp2, 2)))
prec_comp2_corr[1, 2] <- -rho_comp2/(tau_comp2_intercept *
tau_comp2_slope * (1 - pow(rho_comp2, 2)))
prec_comp2_corr[2, 1] <- prec_comp2_corr[1, 2]
mu_comp2_intercept ~ dnorm(mu_comp2_intercept_mean, pow(mu_comp2_intercept_sd,
-2))
tau_comp2_intercept ~ dlnorm(tau_comp2_intercept_meanlog,
pow(tau_comp2_intercept_sdlog, -2))
mu_comp2_slope ~ dnorm(mu_comp2_slope_mean, pow(mu_comp2_slope_sd,
-2))
tau_comp2_slope ~ dlnorm(tau_comp2_slope_meanlog, pow(tau_comp2_slope_sdlog,
-2))
}There are still some minor differences:
decider, the \(\eta\) parameter is part of the 5-parameter
hierarchical vector and has its own hypermean and heterogeneity. In
crmPack, the \(\eta\)
parameter is not part of the hierarchical vector. When there is only
combo arm then this does not matter.Apart from these, the probabilistic models are equivalent. We could also see that in the fit results which are very close to each other.