From 329993b64a871e473f19e411172e07559ebd0498 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 4 Apr 2024 22:08:07 +0200 Subject: [PATCH] 5: Add introduction vignette (#20) --- docs/Project.toml | 6 ++ docs/make.jl | 26 ++++++- docs/src/assets/citations.css | 28 ++++++++ docs/src/index.md | 21 +++++- docs/src/introduction.md | 109 +++++++++++++++++++++++++++++ docs/src/refs.bib | 10 +++ docs/src/small_current_dataset.csv | 101 ++++++++++++++++++++++++++ test/large_historic_dataset.csv | 2 +- test/small_historic_dataset.csv | 2 +- test/test_meta_analysis.jl | 4 +- 10 files changed, 301 insertions(+), 8 deletions(-) create mode 100644 docs/src/assets/citations.css create mode 100644 docs/src/introduction.md create mode 100644 docs/src/refs.bib create mode 100644 docs/src/small_current_dataset.csv diff --git a/docs/Project.toml b/docs/Project.toml index dfa65cd..7be67d5 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,2 +1,8 @@ [deps] +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +SafetySignalDetection = "ad9a3248-762c-4ade-a843-8f185e2c18ff" diff --git a/docs/make.jl b/docs/make.jl index c8f29ae..9a199d7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,15 +1,35 @@ using Documenter using SafetySignalDetection +using DocumenterCitations + +bib = CitationBibliography( + joinpath(@__DIR__, "src", "refs.bib"); + style = :authoryear +) + +module_dir = pkgdir(SafetySignalDetection) +gh = Remotes.GitHub("openpharma", "SafetySignalDetection.jl") +remotes = Dict(module_dir => (gh, "ff9de39")) makedocs( sitename = "SafetySignalDetection", - format = Documenter.HTML(), - modules = [SafetySignalDetection] + format = Documenter.HTML( + assets = String["assets/citations.css"], + size_threshold = nothing + ), + modules = [SafetySignalDetection], + pages = [ + "Home" => "index.md", + "Introduction" => "introduction.md" + ]; + plugins = [bib], + remotes = remotes ) # Documenter can also automatically deploy documentation to gh-pages. # See "Hosting Documentation" and deploydocs() in the Documenter manual # for more information. deploydocs( - repo = "https://github.com/openpharma/SafetySignalDetection.jl.git" + repo = "https://github.com/openpharma/SafetySignalDetection.jl.git", + push_preview = true ) diff --git a/docs/src/assets/citations.css b/docs/src/assets/citations.css new file mode 100644 index 0000000..2d2be76 --- /dev/null +++ b/docs/src/assets/citations.css @@ -0,0 +1,28 @@ +.citation dl { + display: grid; + grid-template-columns: max-content auto; +} + +.citation dt { + grid-column-start: 1; +} + +.citation dd { + grid-column-start: 2; + margin-bottom: 0.75em; +} + +.citation ul { + padding: 0 0 2.25em 0; + margin: 0; + list-style: none !important; +} + +.citation ul li { + text-indent: -2.25em; + margin: 0.33em 0.5em 0.5em 2.25em; +} + +.citation ol li { + padding-left: 0.75em; +} \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 0fbde96..31e0ac8 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,7 +1,26 @@ # SafetySignalDetection.jl -Documentation for SafetySignalDetection.jl +This package implements Bayesian safety signal detection as proposed by [Brock:2023](@citet) +using the [Turing.jl](https://github.com/TuringLang/Turing.jl) framework. + +## Installation + +Install this package with: + +```julia +using Pkg +Pkg.add("SafetySignalDetection") +``` + +## Getting started + +Please have a look at the [Introduction](@ref) vignette to get started. + +## API ```@autodocs Modules = [SafetySignalDetection] ``` + +```@bibliography +``` \ No newline at end of file diff --git a/docs/src/introduction.md b/docs/src/introduction.md new file mode 100644 index 0000000..2e845f1 --- /dev/null +++ b/docs/src/introduction.md @@ -0,0 +1,109 @@ +# Introduction + +We will introduce the use of the `SafetySignalDetection.jl` package with a small example. + +## Data example + +We use a small example data set. + +Note that: + +- The outcome `y` needs to be a bool vector. A different type of vector (e.g. integer) will fail later for the model fitting. +- The column `time` needs to be a `Float64` (time until adverse event or until last treatment or follow up). +- The column `trialindex` needs to be `Int64` (index of historical trials, starting from 1 and consecutively numbered). + +```@example 1 +using CSV +using DataFrames +using SafetySignalDetection + +module_dir = pkgdir(SafetySignalDetection) +csv_path = joinpath(module_dir, "test", "small_historic_dataset.csv") +df = DataFrame(CSV.File(csv_path)) +df.y = map(x -> x == 1, df.y) +``` + +## Meta-analytic prior sampling + +For both parameters $a$ and $b$, we can use beta distributions, specified via `Distributions.Beta()`. + +For the mean $a$, we can set the prior mean according to the background rate of the adverse event. Note that this is per unit of time. So depending how you receive the incidence rate in which units, you need to convert this back to the time unit you are using for the analysis. +The overall sum of the two beta distribution's parameters can be kept constant at a low number. + +For the discount factor $b$, we can set the prior for it closer to 1 for more homogenous trials, +and closer to 0 for more heterogenous trials. +For both priors, the prior information should be low, i.e. the sum of the two parameters +of the beta distribution should be not high, e.g. around 10. + +Note: + +- The number of samples (e.g. 10,000) and the number of threads (1) need to be passed as integers. Otherwise the call will fail. +- We don't need to discard a burn in period now because the NUTS sampler is used internally and is discarding a burn in automatically already. + +```@example 1 +using Distributions + +prior_a = Beta(1 / 3, 1 / 3) +prior_b = Beta(5, 5) + +map_samples = meta_analytic_samples(df, prior_a, prior_b, 10_000, 1) +map_samples[1:5] +``` + +The resulting $\pi^{*}$ samples in `map_samples` are from the meta-analytic prior for the adverse event rate per unit of time in the control arm in the ongoing blinded trial. + +## Closed form approximation + +Now we give the $\pi^{*}$ samples to the `fit_beta_mixture()` function, together with the +number of components. Usually 2 components will be sufficient to achieve a good approximation. +Internally, the classic Expectation Maximization (EM) algorithm is used. + +We can check the approximation graphically by overlaying the resulting probability density function with the +samples histogram. + +```@example 1 +using Plots + +map_approx = fit_beta_mixture(map_samples, 2) + +stephist(map_samples, label = "prior samples", norm = :pdf) +xvals = range(minimum(map_samples), maximum(map_samples); length = 100) +plot!(xvals, pdf(map_approx, xvals), label = "approximate prior") +``` + +## Blinded trial analysis + +Now we can proceed to analyzing the blinded clinical trial. + +The important new ingredients here are: + +- The prior on the experimental adverse event rate. This can be uninformative, or weakly informative using a diluted background rate distribution. +- The proportion of patients in the experimental treatment arm. + +```@example 1 +prior_exp = Beta(1, 1) +prior_ctrl = map_approx +exp_proportion = 0.5 + +csv_current_path = joinpath(module_dir, "docs", "src", "small_current_dataset.csv") +df_current = DataFrame(CSV.File(csv_current_path)) +df_current.y = map(x -> x == 1, df_current.y) + +result = blinded_analysis_samples(df_current, prior_exp, prior_ctrl, exp_proportion, 10_000, 1) +first(result, 5) +``` + +## Summary + +We can now e.g. look at the posterior probability that the rate in the experimental arm, $\pi_E$, is larger than the event rate in the control arm, $\pi_C$. + +```@example 1 +mean(result[!, :pi_exp] .> result[!, :pi_ctrl]) +``` + +We can also create corresponding plots. + +```@example 1 +stephist(result[!, :pi_exp], label = "experimental arm", norm = :pdf) +stephist!(result[!, :pi_ctrl], label = "control arm", norm = :pdf) +``` diff --git a/docs/src/refs.bib b/docs/src/refs.bib new file mode 100644 index 0000000..ee9a8c4 --- /dev/null +++ b/docs/src/refs.bib @@ -0,0 +1,10 @@ +@article{Brock:2023, + author = {Brock, Kristian and Chen, Chen and Ho, Shuyen and Fuller, Greg and Woolfolk, Jared and Mcshea, Cindy and Penard, Nils}, + title = {A Bayesian method for safety signal detection in ongoing blinded randomised controlled trials}, + journal = {Pharmaceutical Statistics}, + volume = {22}, + number = {2}, + pages = {378-395}, + doi = {10.1002/pst.2278}, + year = {2023} +} diff --git a/docs/src/small_current_dataset.csv b/docs/src/small_current_dataset.csv new file mode 100644 index 0000000..8acba6f --- /dev/null +++ b/docs/src/small_current_dataset.csv @@ -0,0 +1,101 @@ +"","y","time","tmt" +"1",0,0.508060037142824,2 +"2",0,0.5208363096442,2 +"3",0,0.556331578741567,2 +"4",0,0.563791380921675,1 +"5",0,0.571754800768225,1 +"6",0,0.590052150745623,1 +"7",0,0.591040918682293,2 +"8",0,0.595221235355885,2 +"9",0,0.597784269735739,1 +"10",0,0.607071721277385,1 +"11",0,0.611330809521675,2 +"12",1,0.613947846980951,2 +"13",0,0.618876835716143,2 +"14",1,0.640354126490284,2 +"15",0,0.645504716144507,2 +"16",0,0.648776558221881,2 +"17",0,0.653263780100834,2 +"18",0,0.664869826653364,2 +"19",0,0.671171271026457,2 +"20",0,0.674994900783847,1 +"21",1,0.679323147659513,2 +"22",0,0.681681905904048,2 +"23",0,0.693983536523541,1 +"24",0,0.713401719711471,1 +"25",0,0.715330274447549,1 +"26",1,0.715554267597332,2 +"27",0,0.724320719842258,1 +"28",0,0.743751407473829,2 +"29",0,0.744621157756557,2 +"30",0,0.749986377628643,2 +"31",0,0.756438894458238,2 +"32",1,0.769038580022455,1 +"33",0,0.781380031033633,2 +"34",0,0.782826549365607,2 +"35",0,0.78411009277818,2 +"36",1,0.791492904824573,2 +"37",1,0.801637802299991,1 +"38",0,0.811791839104605,2 +"39",1,0.837362781844232,1 +"40",0,0.841242652393076,1 +"41",0,0.862928944447668,2 +"42",0,0.869636503880502,1 +"43",0,0.886621511367525,2 +"44",0,0.889937450386028,2 +"45",0,0.89081983407345,2 +"46",0,0.894065672679249,2 +"47",0,0.895545855571809,1 +"48",0,0.909240301239177,2 +"49",0,0.919524707524235,2 +"50",0,0.924650129098574,2 +"51",0,0.945304040451132,2 +"52",0,0.946978089384413,1 +"53",0,0.951859002214535,1 +"54",0,0.957347076003118,2 +"55",0,0.960843806894175,1 +"56",1,0.965779375707409,2 +"57",0,0.966733577783565,1 +"58",0,0.96968774609861,2 +"59",0,0.972657807600261,2 +"60",1,0.979699684335007,1 +"61",0,0.979834431354125,2 +"62",0,0.990429732491942,2 +"63",0,0.993043940746905,2 +"64",0,0.9954099067624,2 +"65",0,0.997159593415083,2 +"66",0,1.00392896736284,1 +"67",0,1.01355110456085,2 +"68",0,1.04904199831044,2 +"69",0,1.06484850264858,1 +"70",0,1.06794773782012,2 +"71",1,1.0712293948771,2 +"72",0,1.08236479130322,1 +"73",0,1.0874704878659,2 +"74",0,1.08815656886374,1 +"75",0,1.09446561418411,2 +"76",1,1.11639969562923,2 +"77",1,1.13149712565995,1 +"78",1,1.15531456740156,2 +"79",0,1.17462593214085,2 +"80",1,1.18276276532017,2 +"81",0,1.18776627686728,1 +"82",0,1.20208448878763,2 +"83",1,1.21019494540132,2 +"84",0,1.21124169832635,1 +"85",0,1.2113759552136,1 +"86",0,1.21897941682307,1 +"87",0,1.22070786681776,2 +"88",0,1.24966304313974,1 +"89",0,1.2527988999829,1 +"90",0,1.25518082340978,2 +"91",0,1.25891242750253,2 +"92",0,1.26357294799558,2 +"93",0,1.26469715973161,1 +"94",0,1.27593270482697,2 +"95",0,1.28273087209758,2 +"96",0,1.29175769180613,1 +"97",0,1.31240415124692,2 +"98",1,1.31259288685395,1 +"99",0,1.33657055292152,2 +"100",1,1.34760702875274,2 diff --git a/test/large_historic_dataset.csv b/test/large_historic_dataset.csv index 9f65459..428483f 100644 --- a/test/large_historic_dataset.csv +++ b/test/large_historic_dataset.csv @@ -1,4 +1,4 @@ -"","y","trial","time" +"","y","trialindex","time" "1",0,1,0.508060037142824 "2",0,7,0.5208363096442 "3",0,2,0.556331578741567 diff --git a/test/small_historic_dataset.csv b/test/small_historic_dataset.csv index bf51bee..7517dbc 100644 --- a/test/small_historic_dataset.csv +++ b/test/small_historic_dataset.csv @@ -1,4 +1,4 @@ -"","y","trial","time" +"","y","trialindex","time" "1",0,5,0.508060037142824 "2",0,4,0.5208363096442 "3",0,2,0.556331578741567 diff --git a/test/test_meta_analysis.jl b/test/test_meta_analysis.jl index 11d55a7..56d0c53 100644 --- a/test/test_meta_analysis.jl +++ b/test/test_meta_analysis.jl @@ -32,7 +32,7 @@ end rng = StableRNG(123) map_small = sample( rng, - meta_analysis_model(df_small.y, df_small.time, df_small.trial, + meta_analysis_model(df_small.y, df_small.time, df_small.trialindex, prior_a, prior_b), NUTS(0.65), 10_000 @@ -54,7 +54,7 @@ end rng = StableRNG(123) map_large = sample( rng, - meta_analysis_model(df_large.y, df_large.time, df_large.trial, + meta_analysis_model(df_large.y, df_large.time, df_large.trialindex, prior_a, prior_b), NUTS(0.65), 10_000