Bayesian Modeling: Yes, No and Maybe Responses

Question

Respondents replied in the following way:

Yes: they will be attending

No: they won't be attending

Maybe: they attach a percentage certainty as an estimate that they'll be attending. E.g. 40% sure I'll make it or 65% sure I'll make it.

I've recently read some literature on Bayesian modeling. Far from knowledgeable in this field but keen to use it for this example and learn in the process.

I intend for the model to be pessimistic so it can account for the worst as the goal of the model is to give an idea of how many people will be attending.

1. The count data observed is modeled by a binomial distribution. I choose it because it is discrete and accounts for there only being two observable states : attends and doesn't attend. So that leaves us with parameters N and p.
2. Modeling p is where it gets difficult. We model p in 3 ways.
3. To conserve the pessimism, we assume that even though people say they are coming there is still between 50% to 100% chance that they will attend. Therefore $p_{1} \sim$ Uniform($0.5,1$).
4. The idea behind the percentage given for maybe-repsonse is to use the trial roulette method of eliciting expert priors. So we end up with some data to which we can fit a normal distribution to. Hence $p_{2} \sim$ Norm($\mu,\sigma$).
5. Cool. I think I'm making sense....perhaps. Now we assume the worst of people such that if a response of no, then they definitely will not attend with 100% probability. Should I model this just as a deterministic value: $p_{3} = 0$

How would I go about putting all this together? I'm using pymc to try and model this. Should I create my own distribution function for $p$? E.g.

$$ p(c) = \begin{cases} U(0.5,1) & \quad \text{if } c_{i} = 1 \\\mathcal{N}(\mu,\sigma) & \quad \text{if } 0<c_{i}<1 \\ 0 & \quad \text{if }c_{i} = 0 \end{cases} $$ Where $c_{i}$ is the set of responses and $c_{i} = 1$ being a yes, 0 a No and and in between value a Maybe.

Would you tell me if what I am doing is absolute nonsense and if I'm on the right path would you please guide me in the right direction. Thank you for your time.

Answer 1

First things first. You basically have two variables, a binary $Y$ that records if they attended ($Y=1$), or not ($Y=0$) and a continuous $X \in [0, 1]$ where the closer $X$ is to $0$, the more certain they are about not attending, and the closer to $1$, the more certain they are that they'd attend. You correctly identified the binomial distribution parametrized by number of possible attendees $n$ and probability of attending $p$, as a likelihood function. In your case $n$ is known, since it is your sample size and $p$ is to be estimated from the data. The $p$ parameter depends on $X$, what leads us to a logistic regression model. First possible choice of a model for your data is simple model assuming continuous $X$:

$$ \begin{align} y_i &\sim \mathsf{Bin}(n, p_i) \\ p_i &= \mathrm{logit}^{-1}\big(\beta_0 + \beta_1 x_i\big) \end{align} $$

Alternatively you may assume something closer to your original description, where we would use a dummy variables: $I_{X=1}$ that returns $1$ if $X=1$ and $0$ otherwise, and $I_{X\in (0,1)}$ that returns $1$ if $X$ is in the $(0, 1)$ range (i.e. excluding the exact $0$ and $1$) and $0$ otherwise. The model would be:

$$ \begin{align} y_i &\sim \mathsf{Bin}(n, p_i) \\ p_i &= \mathrm{logit}^{-1}\big(\beta_0 + \beta_1 I_{X=1} + \beta_2 I_{X\in (0,1)} x_i\big) \end{align} $$

i.e.

$$ \mathrm{logit}(p_i) = \begin{cases} \beta_0 & X = 0 \\ \beta_0 + \beta_1 & X = 1 \\ \beta_0 + \beta_2 x_i & X \in (0, 1) \end{cases} $$

Next, you need to assume priors for $\beta_0,\beta_1,\beta_2$. Assuming that $\Pr(p_i = 0 \mid X = 0) = 0$ does not sound like a reasonable choice since it will lead for your model to always conclude that probability of attending is zero if someone said he will not attend, since $z \times 0 = 0$ no matter what $z$ is. Such prior would it make impossible for your model to learn anything from the data where $X=0$. Same with making strong assumptions like something to leads that $p_i \ge 0.5$ for $X=1$, since there is no objective reason to assume that if someone said that he will do something, then the probability that he will do is $\ge 0.5$ (how many people do you know that told you that they will stop smoking?).

For $\beta_1,\beta_2$ you seemed to assume that they need to be "probabilities". They do not in the model I described, since they can be anything on the real line. Defining such constrains for logistic regression model would make it overtly complicated. In your description you seem to be assuming that the priors are on the $p_i$ values, but this is incorrect since you are interested in the $p_i$ that depends on $x_i$, so you need to include information about such dependence in your model and logistic regression is the most direct and commonly used approach for such cases.

Since logistic function transforms each $z < 0$ to $0 < p < 0.5$ and $z > 0$ to $0.5 < p < 1$, then if you want to choose informative priors for $\beta_0,\beta_1,\beta_2$, use such distributions that lean the resulting values towards positive or negative outcomes, depending on your assumptions.

Answer 2

I think there is something to your approach. However there is a problem with allowing probabilities to be generated by a normal distribution because that can produce values less than zero or greater than one. So I would modify the way in which the probability of attendance depends on the RSVP response. Before doing so, I would like to modify your notation a bit to align it better with what is commonly used in Bayesian analysis.

Let $y = (y_1, \ldots, y_n)$ denote the dataset where $y_i \in [0,1]$ denotes the RSVP response of individual $i$. Let $z_i = 1$ if individual $i$ actually attends and $z_i = 0$ otherwise. Let $x = \sum_{i=1}^n z_i$. I presume you interested in $p(x|y)$, which is the distribution of total attendance given what you observe. (I use "$p$" to denote probability mass or probability density as appropriate.)

Let $\theta_i$ denote the probability that individual $i$ actually attends. Thus \begin{equation} p(z_i|\theta_i) = \textsf{Bernoulli}(z_i|\theta_i) = \begin{cases} \theta_i & z_i = 1 \\ 1 - \theta_i & z_i = 0 \end{cases}. \end{equation}

Now we specify the way in which $\theta_i$ depends on $y_i$. Here is where I modify what you propose. I suggest letting \begin{equation} p(\theta_i|y_i) = \textsf{Uniform}(\theta_i|y_i/2,y_i) , \end{equation} where $\textsf{Uniform}(\theta_i|0,0)$ should be interpreted as a point mass located at $\theta_i = 0$. Note that this distribution is the same as yours for $y_i = 0$ and $y_i = 1$. Moreover, I think it captures the spirit of your "pessimism" for $y_i \in (0,1)$ as well.

At this point we can streamline the model by integrating out the intermediate variable $\theta_i$. In particular, \begin{equation} p(z_i|y_i) = \int_0^1 p(z_i|\theta_i)\,p(\theta_i|y_i)\,d\theta_i = \textsf{Bernoulli}\Big(z_i\,\Big|\,\frac{3\,y_i}{4}\Big) . \end{equation} Now we see that $x$ is the sum of $n$ random variables, each with a Bernoulli distribution. If all the Bernoulli distributions had a common probability of success then the sum would have a binomial distribution; but they don't and so it doesn't.

Nevertheless it is simple to approximate the distribution for $x$ using simulation. Sophisticated software packages are not required for this problem. (And anyway I don't know pymc.) For each simulation $r = 1, \ldots, R$, let $x^{(r)} = \sum_{i=1}^n z_i^{(r)}$, where $z_i^{(r)}$ is a random draw from $\textsf{Bernoulli}(.75\,y_i)$. The larger you make $R$, the more accurate your approximation will be.