Instrumental Variables model in brms?

Hi all,

I am wondering whether there is a way to implement the following instrumental variable model in brms (potentially using the multivariate models functionality)?

(X_i, Y_i)^T \sim \mathcal{N}_2((\xi_i, \eta_i)^T,\Sigma)


\xi_i = \alpha_0 + \sum_k\alpha_k g_{ik}


\eta_i = \beta_0 + \beta_1 \xi_i

Here \mathcal{N}_2 denotes the bivariate normal with mean and covariance matrix as the first and second argument, respectively. i runs over individuals, X_i is viewed as the “exposure” variable, Y_i as the response and g_{ik} for k=1, ..., K are K instruments.

I guess the crux here is that parameters \alpha_0 and \alpha_k's enter in both components of the mean vector. Thus it is not simply done by using brms with bf(y~1+x)+bf(x~1+g)

If it doesn’t work in this parametrisation, maybe there is a way to reparametrise the model to be able to run it brms (or rstanarm)?

Thank you.

I’m a novice with respect to instrumental variables, but it seems to me this is akin to model b14.6 from this link. The statistical model in that particular section was

\begin{align*} \begin{bmatrix} W_i \\ E_i \end{bmatrix} & \sim \operatorname{MVNormal} \begin{pmatrix} \begin{bmatrix} \mu_{\text W,i} \\ \mu_{\text E,i} \end{bmatrix}, \mathbf S \end{pmatrix} \\ \mu_{\text W,i} & = \alpha_\text W + \beta_\text{EW} E_i \\ \mu_{\text E,i} & = \alpha_\text E + \beta_\text{QE} Q_i \\ \mathbf S & = \begin{bmatrix} \sigma_\text W & 0 \\ 0 & \sigma_\text E \end{bmatrix} \mathbf R \begin{bmatrix} \sigma_\text W & 0 \\ 0 & \sigma_\text E \end{bmatrix} \\ \mathbf R & = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix} \\ \alpha_\text W \text{ and } \alpha_\text E & \sim \operatorname{Normal}(0, 0.2) \\ \beta_\text{EW} \text{ and } \beta_\text{QE} & \sim \operatorname{Normal}(0, 0.5) \\ \sigma_\text W \text{ and } \sigma_\text E & \sim \operatorname{Exponential}(1) \\ \rho & \sim \operatorname{LKJ}(2), \end{align*}

where W_i is the criterion, E_i is the focal predictor, and Q_i is the instrumental variable.

1 Like

Thanks @Solomon, this is very close but I think not exactly equivalent, as I would need

\mu_{W,i} = \alpha_W + \beta_E \mu_{E,i}

instead of

\mu_{W,i} = \alpha_W + \beta_E E_i

In other words, the mean parameter \mu_{W,i} of W_i is linear in the mean parameter \mu_{E,i} of E_i, instead of linear in E_i itself.

Intriguing. Seems like you might be able to do that with the non-linear syntax.

Thanks. Would that mean I have to combine non-linear syntax with multivariate syntax? Can you use non-linear parameter combinations appearing in two different outcomes that are mvbind-ed?

Yes, I think that’s what you’d need to do. Unfortunately, I don’t think I’ve done that before. If it were me, I’d start slow and first try to fit a simple intercepts-only multivariate model with the non-linear syntax and build from there.

Following your suggestion to simplify first, I tried something like that, to start with:

prior1 <- prior(normal(0, 1), nlpar = "alpha0")
IV_brm <- brm(bf(mvbind(x,y)~alpha0, alpha0~1, nl=TRUE), data = df, prior=prior1)

Unfortunately I get the following error

Error: The following priors do not correspond to any model parameter: b_alpha0 ~ normal(0, 1)

Non linear parameters may not yet be shared across multiple univariate models. Please check get_prior for how prior specification looks here.

1 Like