A friend pointed me to the parametrisation by Stacey (1962) which is rather neat: if \gamma\sim\mathrm{Gamma}(q^{-2},1), then x=\exp\left(\mu+\sigma/q \log(q^2 \gamma)\right) follows a generalised gamma distribution, and there is an easy mapping between the (q,\mu,\sigma) parametrisatoin and the (a, b, c) parameterization (more here). This parameterization also allows a setup akin to the non-centered parametrisatoin discussed in the section on reparametrization in the documentation by using \gamma rather than x as the latent parameter (although q and \gamma of course remain dependent in the hierarchical prior).
With this setup, I can get the hierarchical model to sample somewhat sensibly but still get a number of divergences. I believe these divergences may be a consequence of \sigma and q being poorly identified (cf. the ratio of \sigma/q in the exponential that defines x and the pair plots below). Any input on how to get this model to sample without divergences would be awesome! I tried redefining x=\exp\left(\mu+\sigma \log(q^2 \gamma)\right) in the hope to make q and \sigma more easily identifiable—but no luck. So here’s the whole story…
Context and model description
I have biomarker measurements y_{ij} of specimen j for patient i. The biomarker concentration can be either a positive number (measured in copies per mL) or a negative test, which means that the concentration was below the level of quantification (LOQ) of the assay. So I model the actual concentration as x_{ij} with a left-censored observation model that gives rise to y_{ij} (the level of quantification is known). The distribution of samples from the patient is informed by parameters q, \mu_i, \sigma at the level of each patient i (the parameters defining the shape of the distribution are pooled amongst all patients). The expected concentration \nu_i for patient i is drawn from a population-level generalised gamma distribution with parameters Q, \lambda, \Sigma, and I determine \mu_i by demanding that \langle y_{ij}\rangle=\nu_i for all j.
In more formal notation,
\begin{align}
\nu_i&\sim\mathrm{GeneralisedGamma}(Q,\lambda,\Sigma)\\
\mu_i&=\log\nu_i + \log\Gamma(q^{-2})-\log\Gamma(q^{-2}(1 + q\sigma)) - 2 \sigma\log(q)/q\\
y_{ij}&\sim \mathrm{GeneralisedGamma}(q,\mu_i,\sigma)\\
x_{ij}&=\begin{cases}
y_{ij}&\text{if }y_{ij}\geq\theta_{ij}\\
\emptyset&\text{if }y_{ij}<\theta_{ij}
\end{cases},
\end{align}
where \theta_{ij} is the level of quantification for the corresponding sample, \emptyset denotes a negative test, and the second line guarantees that \langle y_{ij}\rangle=\nu_i. Here’s the model with sort-of-non-centred parametrisation.
functions {
// Log pdf for the generalised gamma distribution according to Stacey (1962).
real gengamma_lpdf(real x, real q, real mu, real sigma) {
real q2 = q ^ 2;
return - (sigma * x ^ (q / sigma) * exp(-mu * q / sigma) + (2 - q2) * sigma * log(q) +
q * (mu - log(x) + q * sigma * (log(sigma * x) + lgamma(1 / q2)))) / (sigma * q2);
}
// Log cdf for the generalised gamma distribution according to Stacey (1962).
real gengamma_lcdf(real x, real q, real mu, real sigma) {
real log_arg = q / sigma * log(x) - mu * q / sigma - 2 * log(q);
return log(gamma_p(1 / q ^ 2, exp(log_arg)));
}
}
data {
// Information about samples and associations between samples and patients.
int<lower=0> num_patients;
int<lower=0> num_samples;
// Lookup table for patients: i = idx[j] is the patient index for sample j.
int idx[num_samples];
// Information about number of positives and negatives for each patient.
int<lower=0> num_samples_by_patient[num_patients];
int<lower=0> num_positives_by_patient[num_patients];
int<lower=0> num_negatives_by_patient[num_patients];
// Concentration measurements and levels of quantification.
vector[num_samples] load;
vector[num_samples] loq;
}
parameters {
// Population level parameters.
real population_loc;
real<lower=0> population_shape;
real<lower=0> population_scale;
// Individual level parameters
real<lower=0> patient_scale;
real<lower=0> patient_shape;
// Random variable for non-centred setup.
vector<lower=0>[num_patients] patient_gamma_;
}
transformed parameters {
// Evaluate the patient mean in terms of the gamma random variable for an almost
// non-centred parametrisatoin.
vector<lower=0>[num_patients] patient_mean = exp(population_loc + population_scale *
log(population_shape ^ 2 * patient_gamma_) / population_shape);
// Contribution to the patient_loc, evaluated once for efficiency.
real loc_contrib_ = lgamma(1 / patient_shape ^ 2) -
lgamma((1 + patient_shape * patient_scale) / patient_shape ^ 2) -
2 * patient_scale * log(patient_shape) / patient_shape;
// Vector to hold the contributions to the target density for each sample.
vector[num_samples] sample_contrib_;
// Evaluate contributions to the target.
for (j in 1:num_samples) {
// Evaluation the location parameter for this sample.
real loc = log(patient_mean[idx[j]]) + loc_contrib_;
if (load[j] > loq[j]) {
// Account for quantitative observations of the concentration in a sample.
sample_contrib_[j] = gengamma_lpdf(load[j] | patient_shape, loc, patient_scale);
} else {
// Handle left-censoring if the concentration is below the level of quantification.
sample_contrib_[j] = gengamma_lcdf(loq[j] | patient_shape, loc, patient_scale);
}
}
}
model {
// Sample the latent gamma parameters that induce the generalised gamma prior for the
// patient mean.
target += gamma_lpdf(patient_gamma_ | 1 / population_shape ^ 2, 1);
// Add the contributions from all the samples.
target += sum(sample_contrib_);
}
On simulated data (gist), the model is able to recover roughly the right input parameters (I haven’t done a proper credible interval cover analysis in light of the divergences). There are still quite a few divergences, especially near the boundaries, and I probably need to put some time into defining sensible priors to tame the divergences. But I was hoping to get the model to sample so I don’t accidentally hide the pathologies by slapping on priors. Anyway, here are the pairs plots with divergences highlighted.
It looks like the curvature variation at the “ends of the population_shape vs population_scale banana” may be the culprit, but I seem to have reached the limits of my reparametrisation skills. Posterior predictive checks seem fine. Any comments would be greatly appreciated!
