# Binomial measurement error model

**URL:** <https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003>\
**Category:** Modeling\
**Tags:** specification, brms\
**Created:** [July 25, 2024, 2:13pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003 "2024-07-25T14:13:28Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [July 25, 2024, 2:13pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/1 "2024-07-25T14:13:28Z")

</div>

I would like to fit a binomial regression model, where both there is measurement error for both the numerator and denominator. The back-story is that the data comes from a survey that has been adjusted (using sampling weights) to produce the counts and denominators. I do not have access to the original data, just the (adjusted) counts and denominators, and the margin of error for each one. The MOE can be converted to a standard error.

My first question is a Stan model specification question. The Stan User’s Guide gives an [example](https://mc-stan.org/docs/2_18/stan-users-guide/bayesian-measurement-error-model.html) of fitting a Bayesian measurement error model, but the key difference is that, since this is a binomial model, the outcome and denominator counts have to be integers. I tried specifying that model as follows:

```no-highlight
data {
  ...
  vector[N] Y_sd;
  vector[N] trials_sd;
}
parameters {
  ...
  int<lower=0> Y_true[N];
  int<lower=0> trials_true[N];
}
model {
  ...
  Y_obs ~ normal(Y_true, Y_sd);
  trials_obs ~ normal(trials_true, trials_sd);
}

```

There are two (potential) problems with this specification:

1. I got an error message: “(Transformed) Parameters cannot be integers.” How do I declare the “true” values of Y and trials and force them to be integers if I can’t declare them as integers?
2. Even if I’m able to declare them as integers, am I allowed to specify the measurement error model as normal when they are integers? I’m guessing not. What would be a more appropriate distribution here? Keep in mind that I really have no idea about the details of how the adjusted counts and denominators were generated, I was just provided with them and the MOE (which I was told was the width of a 90% CI).

My other question is if there is any way to express this model in brms, but I think that might be a little ambitious.

---

<div class="post-metadata">

**Author:** ![mippeqf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mippeqf/32/18539_2.png) [@mippeqf](https://discourse.mc-stan.org/u/mippeqf)\
**Post date:** [July 25, 2024, 3:14pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/2 "2024-07-25T15:14:28Z")

</div>

This does not seem particularly difficult to model, but can you clarify how your data looks like? Your description indicates that you only have an (adjusted) aggregate ratio of a (binary?) response variable, whereas the `[N]` suffix on your data seems to denote individual survey responses. The best way to go about this in the latter case would be to model response probabilities, whereas for the first scenario, I do not understand the modelling task, as that would more or less be a single data point.

---

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [July 25, 2024, 3:24pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/3 "2024-07-25T15:24:28Z")

</div>

Sorry for not being clear, and thanks for looking at this problem.

I do not have access to the original survey data. What I have is a data frame with one row for every observation (a demographic cell), and the following variables:

- E\_num (the estimated numerator)
- SE\_num (the standard error for that estimate)
- E\_den (the estimated denominator)
- SE\_den (the standard error for that denominator)
- Other variables defining the demographic cell (like age group, gender, county, etc.)

The estimates and SEs were generated by another organization, and are based on applying some sort of sampling weights to the raw survey data (which, again, we don’t have access to). All of the estimates (both numerator and denominator) are integers (presumably rounded).

Thanks!

---

<div class="post-metadata">

**Author:** ![mippeqf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mippeqf/32/18539_2.png) [@mippeqf](https://discourse.mc-stan.org/u/mippeqf)\
**Post date:** [July 25, 2024, 5:30pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/4 "2024-07-25T17:30:57Z")

</div>

You could try to retain the integer constraint by nesting Binomial or Poisson distributions, but I hardly believe that this would yield a tractable model. The following is how I would go about it

```stan
parameters{
  trials_true;
  beta;
}
transformed parameters{  
  theta = inv_logit(X*beta);
}
model{
  trials_obs ~ normal(trials_true, trials_sd);
  Y_obs ~ normal(trials_obs*theta, Y_sd);
}
generated quantities{
  Y_true = trials_obs*theta;
}

```

The key thought here is that instead of modelling

```stan
Y_obs ~ normal(binomial(trials_obs,theta), Y_sd);

```

I’ve used the convergence of the binomial to the normal distribution, given that there are a reasonably large number of individuals in each demographic cell. This will of course yield real estimates of Y\_true and trials\_true. If you intend to embed these values into a larger model structure which requires them to be natural numbers, it would be very much relevant to know more about the framework in order to make appropriate choices in the above.

As a side note, it would be important to normalize Y\_sd and trials\_sd so that the 90% coverage corresponds to a 90% interval in the variance parametrization of the respective normal distribution.

---

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [July 25, 2024, 10:29pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/5 "2024-07-25T22:29:01Z")

</div>

Thanks for the reply.

However, one concern is that we actually do have some pretty small cells, and most of the probabilities are expected to be pretty low (a bunch of the `Y_obs` are zero). Wouldn’t that be exactly the type of situation where the normal approximation doesn’t hold?

Do you have a solution that doesn’t rely on the normal approximation?

---

<div class="post-metadata">

**Author:** ![mippeqf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mippeqf/32/18539_2.png) [@mippeqf](https://discourse.mc-stan.org/u/mippeqf)\
**Post date:** [July 27, 2024, 10:07am UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/6 "2024-07-27T10:07:35Z")

</div>

Hm I don’t know a good solution to that. For small Y, the Binomial is better approximated by a Poisson, which you could try to impose on trials\_obs and Y\_obs. The issue with that is that just like the Binomial, the Poisson has a single parameter, so does not lend itself to parametrizing the mean and variance separately. Hacking yourself around that with sth like Y\_obs ~ (Poisson(lamba) - lambda) \* Y\_sd + lambda doesn’t work either because the Poisson only has support on the positive domain.

Please do post an update if you find a satisfactory solution!

---

<div class="post-metadata">

**Author:** ![ssp3nc3r](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ssp3nc3r/32/120_2.png) [@ssp3nc3r](https://discourse.mc-stan.org/u/ssp3nc3r)\
**Post date:** [July 27, 2024, 7:22pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/7 "2024-07-27T19:22:32Z")

</div>

Since you want a binomial likelihood and using standard errors, would modeling the observations as binomial, and the latent parameters within that as theta work for you?

```no-highlight
data {
  int<lower=0> N; // number of observations
  array[N] int<lower=0> Y_obs; // observed success counts
  array[N] int<lower=0> trials_obs; // observed trial counts
  array[N] real<lower=0> Y_sd; // standard errors for Y
  array[N] real<lower=0> trials_sd; // standard errors for trials
}

parameters {
  vector<lower=0>[N] trials_latent; // latent continuous values for trials
  vector<lower=0,upper=trials_latent>[N] Y_latent; // latent continuous values for Y
}

model {
  Y_obs ~ normal(Y_latent, Y_sd);
  trials_obs ~ normal(trials_latent, trials_sd);
  vector[N] theta = Y_latent ./ trials_latent;
  Y_obs ~ binomial(trials_obs, theta );
}

```

---

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [July 29, 2024, 3:36pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/8 "2024-07-29T15:36:18Z")

</div>

@ssp3nc3r Interesting idea, but we also need to incorporate covariates. i.e.,

```no-highlight
inv_logit(theta) = X * beta;

```

I don’t think your specification allows for this?

---

<div class="post-metadata">

**Author:** ![ssp3nc3r](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ssp3nc3r/32/120_2.png) [@ssp3nc3r](https://discourse.mc-stan.org/u/ssp3nc3r)\
**Post date:** [July 29, 2024, 3:48pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/9 "2024-07-29T15:48:13Z")

</div>

Here’s additional structure:

`theta ~ beta_proportion(inv_logit(X * beta), kappa);`

---

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [July 29, 2024, 4:49pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/10 "2024-07-29T16:49:50Z")

</div>

Is this possible when the model already has theta fixed as:

```no-highlight
vector[N] theta = Y_latent ./ trials_latent;

```

?

Also, in general it just seems odd to model the observed data as Binomial, when it’s really the “true” counts that are part of the DGP. Maybe this is the best we’re going to do though.

---

<div class="post-metadata">

**Author:** ![ssp3nc3r](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ssp3nc3r/32/120_2.png) [@ssp3nc3r](https://discourse.mc-stan.org/u/ssp3nc3r)\
**Post date:** [July 29, 2024, 4:54pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/11 "2024-07-29T16:54:15Z")

</div>

In my example, theta isn’t fixed, it’s elements are parameters.

You have several options for modeling, but first you need to know what parameters you are interested in.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 30, 2024, 7:23pm UTC](https://discourse.mc-stan.org/t/binomial-measurement-error-model/36003/12 "2024-07-30T19:23:17Z")

</div>

> [@jgellar](#):
>
> but the key difference is that, since this is a binomial model, the outcome and denominator counts have to be integers.

> [@jgellar](#):
>
> How do I declare the “true” values of Y and trials and force them to be integers if I can’t declare them as integers?

This isn’t possible in Stan. You can’t literally translate the continuous model to the discrete setting because Stan doesn’t support discrete parameters.

You could

1. do what @mippeqf suggests and use continuous approximations.
2. Alternatively, if the max values of the binomial are small, you can marginalize out the latent true values.
3. Use an alternative software package that allows blocking of discrete parameters, like NumPyro or PyMC in Python or NIMBlE in R. (At least I think these would all support the direct translation of the missing data model.)