using Random, Distributions, DataFrames, Turing, StatsPlots, Plots
#using Pkg
#Pkg.add(["DataFrames", "Turing", "StatsPlots", "Plots"])
In the example below, I implement a simple Success vs. Fail one arm trial design, with a flat prior. Trial data follow a binomial distribution. For illustrative purposes, I set a flat prior, Uniform(0, 1), which is also a Beta(1, 1) distribution. It is a conjugate prior for Binomial. Posterior has a known distribution Beta(Num_Success + 1, n - Num_Success + 1). Nonetheless, I will implement Turing.jl numerical estimate of the posterior, and use the closed-form Beta distribution to verify our results.
Proof,
$$P(p|data) = \frac{P(data|p)P(p)}{P(data)}$$$$P(p|data) \propto Binomial(x; n, p)Beta(1, 1)$$$$P(p|data) \propto \binom{n}{k}p^{k}(1-p)^{n-k} \times 1$$$$P(p|data) \propto p^{(k+1)-1}(1-p)^{(n-k+1)-1} \times 1$$$$P(p|data) \sim Beta(k+1, n-k+1) \times 1$$## For a binary endpoint trial design, success or Fail.
## Sampling distribution - Binomial
## Flat Prior distribution - Uniform(0, 1)
## Specify Trial Parameters
# n = number of patients in a single arm
# p = TRUE probability of success rate given this treatment
N = 10
p = 0.6
##
#= This is it. All we have to do here is to define a macro function for Turing.jl
data are passed in as an argument, in this case, x - number of success.
n is passed in as a parameter - number of trials
Inside the function we simply define and write out the distributions as we would by hand.
In this case we define the prior as a m ~ Uniform() distribution
and data x is binomial given the prior
=#
@model function binaryTrial(x::Vector{Int64}, n::Int64)
m ~ Uniform(0, 1)
x ~ Binomial(n, m)
end
#= Because of this intuitive way of defining distributions in Turing.jl, we can go wild and define complex
layers or priors, like so,
A Normal model with unknown mean and variance.
@model function gdemo(x, y)
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
x ~ Normal(m, sqrt(s²))
y ~ Normal(m, sqrt(s²))
end
(I will showcase more complex cases like this in the next tutorial)
=#
## Simulate trial data
# this seed will give data=7 successes
Random.seed!(5)
d = Binomial(N, p)
# generating just one sample, since one trial is the aggregation of multiple Bernoulli -> one binomial
data = rand(d, 1)
## Complete the function by passing the data into it. Now the object "model" contains all model specifications
model = binaryTrial(data, N)
#= In the next step below we just sample the posterior given a specified sampler
Sample four chains using multiple threads, each with 1000 samples.
I'm using No-U-Turn Sampler (NUTS), an extension of Hamilton Monte Carlo (HMC)
Storing results into "chain"
=#
chain = sample(model, NUTS(), 1000)
## Verify that the numerical posterior density matches that of Beta(Num_Success + 1, n - Num_Success + 1)
histogram(chain[:m], bins=30, alpha=0.4, label="A", normalize = true)
plot!(Beta(data[1] + 1, N - data[1] + 1), labels = "Beta($(data[1]+1), $(N - data[1] + 1))")
## Decision Making
# Against a null of current treatment success rate p* = 0.35
# What is the probability that the treat arm has the same or worse success rate Prob(p <= 0.35)?
sum(chain[:m] .<= 0.35) / size(chain[:m], 1)