I use Bayesian random effects models [$y_i \sim bernoulli\_logit(\beta + \alpha_{subj})$ $\alpha_{subj} \sim normal(0, \gamma)$], the $y$ outcome is binary. Part of the subjects have two observations, part of the subjects have only one observation. I experience that when I use ($y_i, subj_i$) data where a large proportion of the subjects have only one observation the estimated standard deviation of the random effect ($\gamma$) increases to the top region of its weakly informative prior. In other words the model starts to explain the observations with the subject characteristics. I'd like to get realistic estimates of $\gamma$ and $beta$.

Is there a way in Stan / PyMC / JAGS (or maybe another platform) to specify that the $\gamma$ should only be estimated from those cases that have more than one observations? Would that way introduce a bias or another problem?

  • $\begingroup$ @Xi'an I think asking general questions and specific questions are both fine, no details I “should provide” exist. But I very much appreciate if you help with a special case of my question. What additional information would you like to see? As I wrote my idea was to estimate $\gamma$ from the subjects with more than one observation, but estimate $\beta$ from all observations. $\endgroup$ – winerd Dec 21 '17 at 16:33

Sure you can. Think of this as of a missing data problem, you have partial data with some missing cases. To estimate the full model you need first to estimate the parameters given the non-missing data cases, then make some predictions for the missing data and use it in your model. This is how algorithms like EM work and it can be incorporated in a Bayesian model. Since the missing data becomes something to be estimated, in Bayesian scenario you can, and will, additionally assume some priors for the unseen values of the data. Writing down the actual Stan code for such case would be tricky, so after you write down the full model, you could ask for help on the Stan users mailing list.

To give you an example, say that you have two variables $X$ and $Y$ and you use a simple regression model $E(Y) = \beta_0 + \beta_1 X$. Your total sample has size $n$ and for $1,...,k$ first cases you have data on $X$, where for the $k+1,...,n$ you don't (let's call them $\tilde X$). So you can think of your model as of consisting of two "parts". First, you'd estimate $\beta_0,\beta_1$ for the submodel $E(Y_{[1,\dots,k]}) = \beta_0 + \beta_1 X_{[1,\dots,k]}$. Second, given the already known values of $\hat \beta_0, \hat \beta_1$, you estimate $\tilde X_{[k+1,\dots,n]}$ in the submodel $E(Y_{[k+1,\dots,n]}) = \hat\beta_0 + \hat\beta_1 \tilde X_{[k+1,\dots,n]}$. Using the EM algorithm, you would do this iteratively, while in the case of Bayesian model you could define a model where $\beta_0,\beta_1,\tilde X_{[k+1,\dots,n]}$ are the "parameters" to be estimated, where you assume some priors for them. The likelihood for the full data would consists of two "parts" for the missing and non-missing cases:

$$ L(\beta_0, \beta_1, \tilde X_{[k+1,\dots,n]}) = \prod_{i=1}^k \mathcal{N}(y_i | \beta_0 + \beta_1 x_i,\; \sigma^2) \times \prod_{i=k+1}^n \mathcal{N}(y_i | \beta_0 + \beta_1 \tilde x_i,\; \sigma^2) $$

Such approach does not have to introduce much "bias" in your model. It makes your model use all the data you have (in the above examples, all $Y$'s and the $1,\dots,k$ $X$'s), rather then using only the non-missing cases. Obviously, you know less about the cases with missing data, then about the rest of your data, so the model would less precisely learn something about this part of your data. Also you need to remember that what you would assume about the missing data would depend on what you know about the non-missing data (imagine extreme case where you have only single non-missing $X$ and $99$ missing cases, then your best guess about the missing cases would be to assume that they are like the single non-missing value, what could lead to strongly biased results). The above example is very simple, but imagine that you have multivariate model, then the data in the non-missing variables would possibly make your guesses about the missing data in the missing-data variables more precise. As you can imagine, $E(Y) = \square + \square \square$ is a worse case in terms of identifiability, then $E(Y) = \square + \square X_1 + \square \square + \square X_3 + \square X_4 $, where "$\square$" stands for something (partially) unknown to be guessed.

  • $\begingroup$ “you have … missing cases” - Which cases do you see as missing in this scenario? Thanks for helping. $\endgroup$ – winerd Dec 20 '17 at 15:06
  • $\begingroup$ @winerd well, those that you don't have the full data on. $\endgroup$ – Tim Dec 20 '17 at 15:06
  • $\begingroup$ I have ($y_i, subj_i$) observations, neither the $y$ outcome, nor the $subj_i$ subject ID is missing for any observation. The problem arises when many subjects have only one observation. $\endgroup$ – winerd Dec 20 '17 at 15:10
  • $\begingroup$ @winerd so for such subjects you are missing the 2nd, 3rd, ... , k-th observations. It doesn't matter if it is missing by design or not. $\endgroup$ – Tim Dec 20 '17 at 15:12
  • $\begingroup$ I can add a ($y_*=missing,subj_i$) observation for all subjects that have only one observation and I can code the model with the categorical partly observed variable in PyMC. But how will that help? I expect the same large $\gamma$, I expect $\alpha_i$s that are huge if $y_i=1$ or tiny if $y_i=0$ and an imputed $y_*$ posterior sample that is mostly equal to $y_i$. Thanks for helping. $\endgroup$ – winerd Dec 20 '17 at 15:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.