Simulated Data

For full reproducibility: every example in this documentation runs on the study generated by docs/src/simulated_data.jl, reproduced below. simulate_crossover_study(; reference, rng) returns (pk_df, diss_df).

20 subjects, crossover design: reference plus three oral test formulations FormA–FormC. FormD has dissolution data only.
Four profiles removed: subject 1 lacks the reference, subjects 2–4 each lack one test formulation.
Reference: IV bolus, dose 100 for every formulation. reference = :ir swaps in an immediate-release oral reference ($K_a$ = 1.5 /h, 1 h lag), used by Oral Reference Design.
Disposition: one compartment, volume 10, elimination rate 0.2 /h, 10 % log-normal between-subject variability on both.
Absorption and dissolution: Weibull, shape 1.5. Dissolution scales 4, 8, 10, 6 h (FormA–FormD); absorption scale 1.2 × the dissolution scale, times a subject factor with 5 % variability. True correlation: $t_\text{vitro} = t_\text{vivo} / 1.2$, fraction absorbed equal to fraction dissolved.
Sampling: nominal grid 0–48 h every 30 min; interior points jittered by up to 18 % of the step and dropped with probability 8 % (at least 88 % kept). Dissolution is sampled the same way.
No measurement noise: values are read off the exact model curves by linear interpolation.
All draws come from the rng argument, a StableRNG(42) by default, so the study is bit-identical on every Julia version and the global RNG is left untouched. Pass another rng for a different realization. The :iv and :ir variants share subjects and sampling schedules.
Each page includes the file from its @setup block; Documenter copies it next to the rendered pages, so the relative include works in the standalone and the integrated build alike.

using InVitroInVivoCorrelation
using DataFrames
using Random
using StableRNGs

function _jitter_with_dropout(
    rng,
    base_times;
    jitter_fraction,
    dropout_probability,
    minimum_keep_fraction,
)
    step = base_times[2] - base_times[1]
    jitter_width = jitter_fraction * step

    obs_times = copy(base_times)
    for j in 2:(length(obs_times) - 1)
        obs_times[j] += (2.0 * rand(rng) - 1.0) * jitter_width
    end

    keep = trues(length(obs_times))
    keep[2:(end - 1)] .= rand(rng, length(obs_times) - 2) .> dropout_probability

    min_keep = max(3, Int(ceil(minimum_keep_fraction * length(obs_times))))
    n_kept = sum(keep)
    if n_kept < min_keep
        dropped = findall(.!keep)
        for idx in randperm(rng, length(dropped))[1:(min_keep - n_kept)]
            keep[dropped[idx]] = true
        end
    end

    return sort(obs_times[keep])
end

# Linear interpolation on the sorted source grid, extrapolating linearly at the ends —
# the same rule the package's curve sources use.
function _sample_on_times(source_times, source_values, obs_times)
    n = length(source_times)
    return map(obs_times) do t
        j = clamp(searchsortedlast(source_times, t), 1, n - 1)
        t0, t1 = source_times[j], source_times[j + 1]
        v0, v1 = source_values[j], source_values[j + 1]
        v0 + (v1 - v0) * (t - t0) / (t1 - t0)
    end
end

function simulate_crossover_study(; reference::Symbol = :iv, rng::AbstractRNG = StableRNG(42))
    reference in (:iv, :ir) || throw(ArgumentError("reference must be :iv or :ir"))

    times = collect(0.0:0.5:48.0)
    n_subjects = 20
    dose = 100.0

    V_pop, k_pop = 10.0, 0.2
    Ka_pop, tlag_pop = 1.5, 1.0

    subj_V = V_pop .* exp.(0.1 .* randn(rng, n_subjects))
    subj_k = k_pop .* exp.(0.1 .* randn(rng, n_subjects))
    subj_abs_var = exp.(0.05 .* randn(rng, n_subjects))

    shape = 1.5
    diss_scales = Dict("FormA" => 4.0, "FormB" => 8.0, "FormC" => 10.0, "FormD" => 6.0)
    abs_scale_means = Dict("FormA" => 4.8, "FormB" => 9.6, "FormC" => 12.0)

    jitter = (; jitter_fraction = 0.18, dropout_probability = 0.08, minimum_keep_fraction = 0.88)

    ref_name = reference === :iv ? "IV" : "IR"

    pk_df = DataFrame(id = Int[], formulation = String[], time = Float64[], conc = Float64[])
    diss_df = DataFrame(formulation = String[], time = Float64[], diss = Float64[])

    for i in 1:n_subjects
        V, k = subj_V[i], subj_k[i]

        ref_conc = if reference === :iv
            (dose / V) .* exp.(-k .* times)
        else
            t_eff = max.(times .- tlag_pop, 0.0)
            (dose / V) .* Ka_pop ./ (Ka_pop .- k) .* (exp.(-k .* t_eff) .- exp.(-Ka_pop .* t_eff))
        end
        ref_times = _jitter_with_dropout(rng, times; jitter...)
        append!(pk_df, DataFrame(
            id = fill(i, length(ref_times)),
            formulation = fill(ref_name, length(ref_times)),
            time = ref_times,
            conc = _sample_on_times(times, ref_conc, ref_times),
        ))

        uir_true = PolyExponential([1.0 / V], [k])
        for form in ("FormA", "FormB", "FormC")
            abs_scale = abs_scale_means[form] * subj_abs_var[i]
            rate = zeros(length(times))
            pos = times .> 0.0
            t_pos = times[pos]
            rate[pos] .= (shape / abs_scale) .* (t_pos ./ abs_scale) .^ (shape - 1.0) .*
                         exp.(-(t_pos ./ abs_scale) .^ shape)
            conc = convolve(times, dose .* rate, uir_true)
            form_times = _jitter_with_dropout(rng, times; jitter...)
            append!(pk_df, DataFrame(
                id = fill(i, length(form_times)),
                formulation = fill(form, length(form_times)),
                time = form_times,
                conc = _sample_on_times(times, conc, form_times),
            ))
        end
    end

    dropped = Set([(1, ref_name), (2, "FormA"), (3, "FormB"), (4, "FormC")])
    filter!(row -> (row.id, row.formulation) ∉ dropped, pk_df)

    for form in ("FormA", "FormB", "FormC", "FormD")
        diss = 1.0 .- exp.(-((times ./ diss_scales[form]) .^ shape))
        diss_times = _jitter_with_dropout(rng, times; jitter...)
        append!(diss_df, DataFrame(
            formulation = fill(form, length(diss_times)),
            time = diss_times,
            diss = _sample_on_times(times, diss, diss_times),
        ))
    end

    return pk_df, diss_df
end