I am attempting to synthesise willingness-to-pay (WTP) estimates from the literature for a given intervention as part of a systematic review. The goals of this analysis are to understand roughly the parameters and moments of the underlying distribution of WTP, as well as how that distribution may differ across countries. The outputs will be used for a few things - a discrete choice experiment, an economic evaluation, and maybe as priors in some future analysis. I’m fairly new to inferential statistics as most of my previous work involves cost-effectiveness research and prediction models, and am <1 year post-PhD, so mea culpa if there is any basic stuff I should know but don’t.
Briefly, the problem is as follows:
There are 8 studies reporting willingness-to-pay using questionnaires in a variety of countries, currencies and years (i.e. both inflation and currency conversion are a concern)
Willingness-to-pay is elicited as an ordinal response variable (e.g., $0, $1 to $100, $100 to $200, etc).
Some studies report an unbounded final value (e.g., “More than $1,000”)
Represented as a factor, there are 28 levels, but the distance between levels varies wildly even within the same study (not much validation going on in willingness-to-pay surveys, apparently).
The approach that seems most intuitive to me is to create an “individual” dataset corresponding to the sample size for each response and upper bound willingness to pay values in each study, convert responses to a common currency/year value, and run a Bayesian ordinal regression with country as a predictor and study ID as a random effect.
I briefly discussed this matter with an economist who suggested I try a Bayesian ordinal cumulative probit using brms. I have some passing familiarity with the package, but have previously used NIMBLE.
My questions:
Is the ordinal cumulative probit reasonable here? Would it be able to natively account for these unbounded responses, or are significant assumptions needed?
If I’ve fully specified the model already, are there some good model diagnostics to run besides trace plots?
How can I report the regression results as distribution parameters?
With 28 levels you are implying that different studies use different WTP categories and you are pooling across studies. I am fairly certain that simple pooling is impossible. Suppose that one study used intervals [50, 100) and [100, 150) and another used [50, 75), [75, 100), [100, 125), [125, 150). You can consider all of the lower/upper value pairs as representing interval-censored values, and get the Turnbull estimate of the cumulative distribution function. This will not have steps at all possible left and right endpoints. If you have covariates this generalizes to a proportional odds model with interval censoring as handled by R rmsb::blrm (Bayesian using Stan) and rms::orm (frequentist).
Whichever approach you take, without having uniformity of intervals there is not reliable way to do this simply.
Thanks Frank, this is more or less what I was afraid of. There are different levels across studies, with some overlap that becomes more complicated when you consider that rarely does $100 mean the same thing between one study and the next due to inflation.
Ultimately I’m interested in just one covariate: the country setting for each study. However, if there is no way to fit them into a single model without torturing the data, I would prefer to just report things descriptively, but still meet the objective of providing a prior for future research.
My last resort is pretty low-rigour: take the number of responses for each level, ascribe them to the upper limit of the interval (assuming roughly equidistant endpoints for the unbounded intervals), and fit a Gamma distribution representing the latent WTP for each study, region, and overall. I’d probably represent this in a ridgeline plot. I’m shying away from this at the moment because there are too many strong assumptions, but I suppose it’s better than nothing.
Why not take the survey results as censored observations of an underlying continuous latent variable, and estimate a model built upon that variable? This way of approaching the problem would have the advantage of demanding that you confront whatever ‘real’ underlying econopsychometric phenomena are at work, because that latent variable brings you conceptually closer to them.
The approach I outlined does that without latent variables, but the latent variable approach (which can also be done with Stan) is a good one too.
I would avoid any method that just uses the upper limits.
A very approximate alternative to all this is to map all the distinct intervals into large intervals and map all the raw data that way. But this loses some information.
Thanks @davidcnorrismd and @f2harrell, this is really helpful. I think the latent continuous variable approach would be my preference, which as noted perhaps better represents the fact that these are approximations of some underlying WTP value.
If you know of any resources that explain how to structure the regression equation to reflect, for example, differences between (100, 200) and the censored 200+ value and how to ensure you’re estimating the underlying latent variable, that would be much appreciated. Happy to use rmsb as well.
I have found JAGS’s ‘observable function’ dinterval models censoring very cleanly. (It’s billed in JAGS manual as a conceptually clarified version of the I(,) construct of BUGS.) My paper linked below shows an application, where the censoring occurs due to categorization of a latent toxicity threshold.
For the non-latent variable approach see this which uses the rms package Ocens function while specifying the dependent variable, e.g. y <- Ocens(c(125, 100, 200), c(125, 200, Inf)) indicates a uncensored 125 value, interval censored 100-200, right censored at 200.
Thank you greatly to both @f2harrell and @davidcnorrismd for your feedback on this. I have a simple model in brms now working with very reasonable findings. In brief, it looks like this:
with a hurdle_gamma() family and 4 chains. The only prior I’ve specified is a half-Cauchy for the intercept, and the rest use the default weakly informative flat priors in brms. The choice of hurdle Gamma is because there is likely a true 0 here where people are ideologically opposed to paying for the treatment.
Trace plots look reasonably good, and there’s a bit of divergence as it explores the long right tail but nothing major. I’m playing around with removing the Setting variable (maybe too much overlap with the random intercepts) and perhaps increasing the thinning to reduce autocorrelation. The posterior predictions look really good and suggest that I have managed to get at the underlying latent willingness-to-pay variable.
I really appreciate the push to do this in a fully Bayesian approach; I’ve never had the chance to do something like this from start to finish, and it’s been quite fun/satisfying. The ability to make valid statements about the probability of an outcome being within a certain range is incredibly useful for decision-making.
When we go to publish, would either of you mind being listed in the acknowledgements?