Hi all, I just started with RxInfer.jl for simulations of patient recruitment in clinical trials. The simplest model is a simple Poisson-Gamma model, where each clinical site draws the (time homogeneous) rate for the Poisson process describing patient recruitment from a Gamma distribution. I am trying to implement in RxInfer but getting ERROR: UndefRefError: access to undefined reference. I am sure this is some simple error on my part but I cannot find it. Appreciate any help.
using Distributions
"""
Follows 3.2.1 in the paper. Simulate data from a clinical trial enrollment process.
Note β is a rate parameter, but Distributions uses scale; we take 1/β in the function.
Note that the paper considers the "time step" to be a day.
"""
function pg_homogeneous(α, β, u, T)
N = length(α)
@assert N == length(β) == length(u)
λ = rand.(Gamma.(α, 1 ./ β)) # shape, rate parameterization
outmat = zeros(Int, N, T)
for i in 1:N
rate = [k < u[i] ? 0 : λ[i] for k in 1:T]
outmat[i, :] = rand.(Poisson.(rate))
end
return cumsum(outmat, dims=2) # return cumulative number of pts
end
# example with time-homogeneous rates
α = 1/(1.2^2)
β = α/0.02
N = 200
T = 400
u = rand(1:29*4, N)
recruitment_target = 1000
out_homogeneous = pg_homogeneous(fill(α, N), fill(β, N), u, T)
final_date = findfirst(sum(out_homogeneous, dims=1) .>= recruitment_target)[2]
out_homogeneous = out_homogeneous[:, 1:final_date]
using RxInfer
# The inference model for the homogeneous PG model.
# This assumes that Gamma(α,β) is the common distribution of the
# rate for all centers (the simulation function in the paper generalizes that somewhat
# but all the math in the paper makes this simplifying assumption).
@model function pg_homogeneous_rxinfer(data, T, N, u)
# these are quite terrible priors for the Gamma(α,β)
# parameters, but the actual conjugate prior is quite a strange non-standard
# distribution (see Gamma entry in https://en.wikipedia.org/wiki/Conjugate_prior#When_likelihood_function_is_a_continuous_distribution)
# so let's use these for testing
α ~ Uniform(0, 10)
β ~ Uniform(1, 150)
# the main stochastic model
for i in 1:N
λ[i] ~ Gamma(shape=α, scale=β) # shape, scale parameterization
for k in 1:T
if k ≥ u[i] # recruitment has started
data[i, k] ~ Poisson(λ[i]) # Poisson draw at each day
end
end
end
end
infer_homogeneous = infer(
model = pg_homogeneous_rxinfer(T=final_date, N=N, u=u),
data = (data = out_homogeneous, )
)
Hey! Sorry I don’t have a definite solution but could this have something do with the conjugacy stuff? I understand that RxInfer requires it by default and if you have a non-conjugate model you need to take some extra steps (?). (Coming from medicine myself, cool application by the way!)
I can’t quite figure out what’s happening. Seems like a malformed Vector value, with uninitialized fields, is constructed somehow, then trying to index into it throws UndefRefError. Either there’s a Julia bug, or (more likely, I guess) a bug in one of the packages here.
Try reporting a bug to RxInfer.jl, perhaps the maintainers figure it out.