4: add blinded analysis model and corresponding sample function (#18)

* add code from daniel branch

* fix with model macro syntax

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")

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)
+
+"""
+@model 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

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,

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")

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

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