Difficulties with a Bayesian formulation of a model for human timing data

Question

The Wing-Kristofferson model is a simple model of the behavior of a human trying to drum out a steady beat (that is, trying to mimic a metronome). Let $y_i$ be the $i$th interval between two drum beats then the model is:

$$y_i \sim \text{Normal}(\mu + m_i - m_{i-1}, \sigma_y) $$ $$m_i \sim \text{Normal}(0, \sigma_m) $$

where $\mu$ is the mean interval of a time keeper and $m_i$ is the error of the $i$th drum stroke. What one wants is to estimate the parameters $\sigma_y$ and $\sigma_m$ ( while $\mu$ is considered more of a nuisance parameter).

I've been trying to implement this model in a Bayesian framework by treating the $m_i$ as missing data. I've got two problems that I can't get my head around and that I would really appreciate some input on:

1. What would be appropriate "vague/objective" priors for the SD parameters $\sigma_y$ and $\sigma_m$?

I've been trying to to find "objective" priors for $\sigma_y$ and $\sigma_m$ that let the data dominate and that are relatively unbiased. Using simulated data I've tried a number of priors, non which have worked well so far. Flat priors on $\sigma_y$ and $\sigma_m$ tend to overestimate $\sigma_m$ but flat priors on $\log(\sigma_y)$ and $\log(\sigma_m)$ tend to pull $\sigma_m$ close to zero.

2. Is there a way to reduce the correlation between the parameters?

As $m_i$ are treated as missing data, and thus are parameters of the model, the number of parameters increase with the number of datapoints $y_i$. This wouldn't be too bad if it wasn't for the fact that both the parameters $m_i$ are heavily correlated both with each other and with $\sigma_m$. Another issue is that that $\sigma_m$ and $\sigma_y$ are heavily correlated. All these correlations makes the models hard to fit using MCMC frameworks such as JAGS or STAN. I've been tinkering around with different reparametizations of the model, but so far I haven't found anything that works. Still when the number of data points are low ($n < 60$) this isn't too much of a problem.

Answer 1

I suggest that instead of treating the $m_i$ as missing data, you integrate them out, since this marginalization leaves the form of the data distribution unchanged (in the sense that it's still a multivariate normal distribution). Marginally of the $m_i$, $\mathrm{E}(y_i)=\mu$ and

$\mathrm{Cov}\left(y_{i},y_{j}\right)=\begin{cases}2\sigma_{m}^{2}+\sigma_{y}^{2}, & i=j,\\-\sigma_{m}^{2}, & \left|i-j\right|=1,\\0, & \mathrm{otherwise.}\end{cases}$

That's a nice sparse covariance matrix right there.

That leaves you with a three-parameter problem, and in fact $\mu$ can be marginalized out too if the prior is conjugate or flat, so now you're down to a 2D posterior distribution. Surely STAN can handle that -- but even if it can't, you can just plot the damn thing, integrate it numerically, and use old reliable.

According to Gelman's Bayesian Data Analysis, 2$^{nd}$ ed., equation 14.14 on page 375, for a flat prior on the mean parameter,

$p(\Sigma_y|y) \propto p(\Sigma_y)|\Sigma_y|^{-1/2}|V_{\beta}|^{1/2}\exp\bigg(-\frac{1}{2}(y-X\hat{\beta})^T\Sigma_y^{-1}(y-X\hat{\beta})\bigg),$

in which

• $y$ is the complete data vector,
• $X$ is, in your case, a vector of ones,
• $\Sigma_y$ is the covariance matrix specified by the above expression for the covariance, (and hence $p(\Sigma_y|y) = p(\sigma^2_m,\sigma^2_y|y)$)

and

$\hat{\beta} = (X^T\Sigma_y^{-1}X)^{-1}X^T\Sigma_y^{-1}y$

$V_{\beta} = (X^T\Sigma_y^{-1}X)^{-1}$

Note that these above expressions make $\hat{\beta}$ and $V_{\beta}$ functions of the unknown parameter $\Sigma_y^{-1}$.

As for the prior $p(\sigma^2_m,\sigma^2_y)$, a good place to start is Gelman's 2006 paper Prior distributions for variance parameters in hierarchical models.