Is appropriate to use empirical Bayes (EB) in this way?

Question

Background. I have data from a study where participants make a series of judgments (a series of decisions with a binomial outcome, either $y=1$ or $y=0$). I have a model of the underlying decision-making process, which has free parameters that can be estimated from the data. My interest is ultimately comparing these parameters across different groups of participants. I have a relatively small number of observations on each participant, so I am using a hierarchical Bayesian approach (I wrote my model in Stan) to pool data across individuals.

The problem. For the comparison to be accurate I need to take into account that some of the responses might be purely "random", in the sense that they are not predictable (e.g. participant is distracted and pushes either one of the two buttons with probability 0.5). In other words, the likelihood of the observed responses is a mixture between complete chance, $\frac{1}{2}$, and the probability predicted by the model, call it $p^*$. The mixture is governed by a parameter $\lambda$ which basically correspond to the probability of giving a "random" response (often referred to as "lapse" of attention). In other words the probability of observing $y=1$ in any trial is given by $\lambda \cdot \frac{1}{2} + \left(1 -\lambda \right) \cdot p^* $.

More formally, the model could be expressed as $$ P\left(y_i=1 \mid j\right) = \frac{\lambda_j}{2} + \left( 1- \lambda_j\right)\cdot\Phi\left[\beta_0 + u_{j0} + \left(\beta_1 + u_{j1} \right)x_i\right] $$ where $j$ indicate the participant, $\beta_0, \beta_1$ are fixed-effects; $u_0, u_1$ are participant-specific random effects; and $\lambda_j$ is the lapse rate of subject $j$.

To have an idea of how frequently these lapses occur, I have included few "catch" trials, that is decisions where the correct answer is so obvious that we can safely assume that participants would make a mistake only if they responded at chance. Thus the frequency of errors in these catch trials could be taken as an estimate of the lapse rate, although I have only have few catch trials (either 6 or 12, depending on the participant). I want to use a multilevel approach also for the parameter $\lambda$, specifically by assuming a beta distribution, $\lambda \sim \text{Beta} \left(a, b \right)$. However, while I can assign sensible hyper-priors to all other parameters, I have difficulty deciding and justifying on an appropriate hyper-prior for $a$ and $b$. Note that I am not interested in the lapse rates per se, I just want to control for the probability of lapses them while estimating and comparing the other parameters across groups.

Questions. I was considering using an approach which I think could be defined as empirical Bayes, where for each participant $j$ I estimate the lapse rate $\hat \lambda_j$ as simply the ratio of the number of errors (in catch trials) to the number of catch trials; then I estimate the values of $a$ and $b$ via MLE (by maximizing the likelihood of $\hat \lambda_j$ under a Beta distribution; see the plot below), and then plugging in these estimates in the multilevel model above, as a prior distribution for the participant-specific lapse rates.

• Is this approach legit?
• I have zero experience in empirical Bayes (EB), but I have the impression that in most cases it is implemented as an iterative process resembling expectation maximization, whereas in my case I wouldn't do any iterations. Is that still a valid way of applying EB? (Pointers to relevant references are appreciated!)

• Is it a problem if in my model I would have some parameters with an "empirical" prior, and some others with standard Bayesian priors and hyper-priors?

  

Answer 1

Whether or not your approach is legitimate depends in large part about how you describe your approach when publishing or presenting your results. If you are completely open about your approach and your process then the reader is able to judge your approach for themselves. I state this because statistics so often involves subjective choices that have no clear right or wrong answer that the best approach is to simply make all of your choices open. What I mean is an approach can be a perfectly legitimate use of Empirical Bayes but the reader might take issue with Empirical Bayes.

For the purposes of your model though: you have picked one approach that is consistent with Empirical Bayes work. See for example here: http://varianceexplained.org/r/empirical_bayes_baseball/

Reference supporting your approach: https://www.jstor.org/stable/2669771?seq=1#page_scan_tab_contents

Again as long as you make it clear to the reader/consumer of your work how you chose each prior distribution you allow the reader to make the decision of how much they agree with your analysis.

Another approach: If this approach is making you uncomfortable then there is an approach that I personally like better and that is to pick a prior distribution that generates reasonable data. This is nicely demonstrated here: https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/rssa.12378 (especially figure 4). Basically what you do is simulate your data from the prior and see if it reasonable approximates the real data. I also recommend reading this: https://arxiv.org/pdf/1708.07487.pdf to get an understanding of the entire thought process.

Answer 2

Let me slightly diverge from your exact question. You are describing your intended model, where the probability of giving the answer $y=1$ is modeled as

$$ p(y_i=1) = \lambda\,0.5 + (1 - \lambda) \, p^*_i $$

Notice that the proposed model can be described in different form

$$ \lambda\;0.5 + (1 - \lambda) \, p^*_i = \alpha + \gamma_i $$

in such case, you could write

$$ \lambda = \frac{\alpha}{0.5}, \qquad p^*_i = \frac{\gamma_i}{1-\lambda} $$

If you think about it, then it seems that you can re-define your model to logistic regression with $\alpha$ and $\gamma$ being unbounded, real-valued parameters, that get turned into probabilities by passing them through the logistic function $\sigma(\cdot)$,

$$ p(y_i=1) = \sigma(\alpha + \gamma\,d_i) $$

where $d$ is an indicator that is equal to $1$ for the regular trials, and $0$ for the "catch" trials. In such case, probability of $y=1$ for the "catch" trial is simply $\sigma(\alpha)$. You should be able to learn $\alpha$ from your data in a single step, without any empirical Bayesian tricks. In such case, $\alpha$ would be the "base rate" for "catch" trials, and the added effect of non-random answers would be modeled by $\gamma$.

Finally, this enables us to re-define your complete model as

$$ p(y_i=1) = \sigma\big(\alpha + d \cdot [\beta_0 + u_{j0} + (\beta_1 + u_{j1}) \,x_i]\big) $$

I do not have any formal arguments to back my thesis, but for me, this formulation seems to be simpler and more flexible (e.g. $\lambda$ and $p^*$ are not constrained to unit interval). You could argue with that, but for me, such formulation is also cleaner in terms of interpretability, because you model the additive effect directly, withouth the weighting by $\lambda$.

As about priors, for the above model you could choose something like $\alpha \sim \mathcal{N}(0, \tau)$, where $\tau$ would control the variability around $0.5$ assumed a priori for $p(y=0)$.