From ce92533f8dec0970af2ed510b542c0d0570e951c Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Mon, 18 Mar 2024 15:23:59 +0100 Subject: [PATCH 1/2] add code from daniel branch --- src/SafetySignalDetection.jl | 3 ++ src/blinded_analysis.jl | 75 +++++++++++++++++++++++++++++++++++ src/meta_analysis.jl | 8 ++-- test/runtests.jl | 1 + test/test_blinded_analysis.jl | 42 ++++++++++++++++++++ test/test_meta_analysis.jl | 2 +- 6 files changed, 126 insertions(+), 5 deletions(-) create mode 100644 src/blinded_analysis.jl create mode 100644 test/test_blinded_analysis.jl diff --git a/src/SafetySignalDetection.jl b/src/SafetySignalDetection.jl index 13b0866..edae280 100644 --- a/src/SafetySignalDetection.jl +++ b/src/SafetySignalDetection.jl @@ -10,10 +10,13 @@ using LinearAlgebra using ExpectationMaximization export + blinded_analysis_model, + blinded_analysis_samples, meta_analysis_model, meta_analytic_samples, fit_beta_mixture +include("blinded_analysis.jl") include("meta_analysis.jl") include("fit_mle.jl") include("fit_beta_mixture.jl") diff --git a/src/blinded_analysis.jl b/src/blinded_analysis.jl new file mode 100644 index 0000000..3bec4f2 --- /dev/null +++ b/src/blinded_analysis.jl @@ -0,0 +1,75 @@ +""" + +Blinded Analysis Model + +This Turing model is used to generate posterior samples of the adverse event probabilities +`pi_exp` in the experimental arm and `pi_ctrl` in the control arm given a blinded analysis +of a trial with `exp_proportion` ratio of experimental arm patients relative to all patients. + +blinded_analysis_model( + y::Vector{Bool}, + time::Vector{Float64}, + prior_exp::Distribution, + prior_ctrl::Distribution, + exp_proportion::Float64) + +""" +function blinded_analysis_model( + y::Vector{Bool}, + time::Vector{Float64}, + prior_exp::Distribution, + prior_ctrl::Distribution, + exp_proportion::Float64) + + n = length(y) + pi_exp ~ prior_exp + pi_ctrl ~ prior_ctrl + mu_exp = log(-log(1 - pi_exp)) + mu_ctrl = log(-log(1 - pi_ctrl)) + + for i in 1:n + x_exp = mu_exp + log(time[i]) + x_ctrl = mu_ctrl + log(time[i]) + prob_exp = 1 - exp(-exp(x_exp)) + prob_ctrl = 1 - exp(-exp(x_ctrl)) + y[i] ~ MixtureModel( + Bernoulli[ + Bernoulli(prob_exp), + Bernoulli(prob_ctrl) + ], + [exp_proportion, 1 - exp_proportion] + ) + end + +end + + +""" + +Blinded Analysis Posterior Samples Generation + +This function wraps the Turing model `blinded_analysis_model` and runs it for a data frame `df` with: + - `y`: `Bool` (did the adverse event occur?) + - `time`: `Float64` (time until adverse event or until last treatment or follow up) + +Note that arguments for the number of samples per chain and the number of chains have to be passed as well. + +It returns a `DataFrame` with the posterior samples for `pi_exp` and `pi_ctrl`. +""" +function blinded_analysis_samples( + df::DataFrame, + prior_exp::Distribution, + prior_ctrl::Distribution, + exp_proportion::Float64, + args... + ) + chain = sample( + blinded_analysis_model(df.y, df.time, prior_exp, prior_ctrl, exp_proportion), + NUTS(0.65), + MCMCThreads(), + args... + ) + df = DataFrame(chain) + select!(df, [:pi_exp, :pi_ctrl]) + +end \ No newline at end of file diff --git a/src/meta_analysis.jl b/src/meta_analysis.jl index dd949dd..477df25 100644 --- a/src/meta_analysis.jl +++ b/src/meta_analysis.jl @@ -43,16 +43,16 @@ meta_analysis_model( y[i] ~ Bernoulli(prob) end -end; +end """ Meta Analytic Prior Samples Generation This function wraps the Turing model `meta_analysis_model` and runs it for a data frame `df` with: - - `y`: Bool (did the adverse event occur?) - - `time`: Float64 (time until adverse event or until last treatment or follow up) - - `trialindex`: Int64 (index of trials, starting from 1 and consecutively numbered) + - `y`: `Bool` (did the adverse event occur?) + - `time`: `Float64` (time until adverse event or until last treatment or follow up) + - `trialindex`: `Int64` (index of trials, starting from 1 and consecutively numbered) meta_analytic_samples( df::DataFrame, diff --git a/test/runtests.jl b/test/runtests.jl index 7f641e0..fb3c5a5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ using Turing using SafetySignalDetection include("test_helpers.jl") +include("test_blinded_analysis.jl") include("test_meta_analysis.jl") include("test_fit_beta_mixture.jl") include("test_fit_mle.jl") diff --git a/test/test_blinded_analysis.jl b/test/test_blinded_analysis.jl new file mode 100644 index 0000000..c744046 --- /dev/null +++ b/test/test_blinded_analysis.jl @@ -0,0 +1,42 @@ +@testset "blinded_analysis_model works as expected" begin + rng = StableRNG(123) + + n_per_arm = 50 + df = DataFrame( + y = vcat( + rand(rng, Bernoulli(0.2), n_per_arm), + rand(rng, Bernoulli(0.5), n_per_arm) + ), + time = rand(rng, Exponential(1), 2 * n_per_arm) + ) + + chain = sample( + rng, + blinded_analysis_model(df.y, df.time, Beta(2, 8), Beta(9, 10), 0.5), + NUTS(0.65), + 1000 + ) + + check_numerical(chain, [:pi_exp], [0.091], rtol=0.001) + check_numerical(chain, [:pi_ctrl], [0.659], rtol=0.001) +end + +@testset "blinded_analysis_samples works as expected" begin + rng = StableRNG(123) + + n_per_arm = 50 + df = DataFrame( + y = vcat( + rand(rng, Bernoulli(0.2), n_per_arm), + rand(rng, Bernoulli(0.5), n_per_arm) + ), + time = rand(rng, Exponential(1), 2 * n_per_arm) + ) + + samples = blinded_analysis_samples(df, Beta(2, 8), Beta(9, 10), 0.5, 100, 1) + + @test typeof(samples) == DataFrame + @test nrow(samples) == 100 + @test ncol(samples) == 2 + @test names(samples) == ["pi_exp", "pi_ctrl"] +end \ No newline at end of file diff --git a/test/test_meta_analysis.jl b/test/test_meta_analysis.jl index 035ced7..11d55a7 100644 --- a/test/test_meta_analysis.jl +++ b/test/test_meta_analysis.jl @@ -1,4 +1,4 @@ -@testset "Check that meta_analysis_model works as before" begin +@testset "meta_analysis_model works as before" begin rng = StableRNG(123) n_trials = 5 From a526a463642b6f315cd54bc1484a8d5a4d25f97f Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Mon, 18 Mar 2024 17:59:32 +0100 Subject: [PATCH 2/2] fix with model macro syntax --- src/blinded_analysis.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blinded_analysis.jl b/src/blinded_analysis.jl index 3bec4f2..407f12b 100644 --- a/src/blinded_analysis.jl +++ b/src/blinded_analysis.jl @@ -14,7 +14,7 @@ blinded_analysis_model( exp_proportion::Float64) """ -function blinded_analysis_model( +@model function blinded_analysis_model( y::Vector{Bool}, time::Vector{Float64}, prior_exp::Distribution,