Partial pooling of coefficients in a temporal DAG model

Interfaces brms
#1

Hi all,

I want to use brms to model a DAG like the one below:

\begin{equation} \begin{matrix} X_1 & \rightarrow & X_2 & \rightarrow & X_3 \\ \downarrow & & \downarrow & & \downarrow \\ Y_1 & \rightarrow & Y_2 & \rightarrow & Y_3 \\ \end{matrix} \end{equation}

Where the subscripts \{1, 2, 3\} represent time.

Suppose this is my data:

print(XY)

     X1    X2    X3    Y1    Y2    Y3
1 108.2 110.9 131.0 162.2 143.1  89.1
2  89.8  88.1  87.7 156.5 248.6 287.6
3  97.7 102.7 125.5 235.0 415.4 657.1
4  98.4  84.1  84.6 226.3 375.5 506.0
5  85.1  97.2 125.3 202.6 274.9 269.6

Where the model looks something like this:

\begin{equation} X_{3i} \sim N(\alpha_2 X_2, 10) \\ X_{2i} \sim N(\alpha_1 X_1, 10) \\ X_1 \sim N(100, 10) \\ Y_{3i} \sim N(\beta_2 Y_{2i} \space + \space \lambda_3 X_{3i}, 45) \\ Y_{2i} \sim N(\beta_1 Y_{1i} \space + \space \lambda_2 X_{2i}, 45) \\ Y_{1i} \sim N(\lambda_1 X_{1i}, 45) \\ [ \alpha_1, \alpha_2 ]^T \sim MVN(0, \Sigma_{\alpha}) \\ [ \beta_1, \beta_2 ]^T \sim MVN(0, \Sigma_{\beta}) \\ [ \lambda_1, \lambda_2, \lambda_3 ]^T \sim MVN(0, \Sigma_{\lambda}) \\ \Sigma_{\alpha} = \begin{bmatrix} \sigma_{\alpha_1} & 0 \\ 0 & \sigma_{\alpha_2} \end{bmatrix} R_{\alpha} \begin{bmatrix} \sigma_{\alpha_1} & 0 \\ 0 & \sigma_{\alpha_2} \end{bmatrix} \\ \Sigma_{\beta} = \begin{bmatrix} \sigma_{\beta_1} & 0 \\ 0 & \sigma_{\beta_2} \end{bmatrix} R_{\beta} \begin{bmatrix} \sigma_{\beta_1} & 0 \\ 0 & \sigma_{\beta_2} \end{bmatrix} \\ \Sigma_{\lambda} = \begin{bmatrix} \sigma_{\lambda_1} & 0 & 0 \\ 0 & \sigma_{\lambda_2} & 0 \\ 0 & 0 & \sigma_{\lambda_3} \end{bmatrix} R_{\lambda} \begin{bmatrix} \sigma_{\lambda_1} & 0 & 0 \\ 0 & \sigma_{\lambda_2} & 0 \\ 0 & 0 & \sigma_{\lambda_3} \end{bmatrix} \\ \sigma_{\alpha_1}, \sigma_{\alpha_2} \sim Exponential(1) \\ \sigma_{\beta_1}, \sigma_{\beta_2} \sim Exponential(1) \\ \sigma_{\lambda_1}, \sigma_{\lambda_2}, \sigma_{\lambda_3} \sim Exponential(1) \\ R_{\alpha} \sim LKJ(1) \\ R_{\beta} \sim LKJ(1) \\ R_{\lambda} \sim LKJ(1) \end{equation}

Basically, the aim is to do partial pooling for the \alpha, \beta, \lambda coefficients across time, and to understand their correlation across time also. What I am struggling to do in brms is passing in multiple dataframes, say, a dataframe each for X_1, X_2, X_3, Y_1, Y_2, Y_3 so that each dataframe can have another column to indicate the time index and the dependent variable for that node.

So far, because I do not know how to pass in multiple dataframes into brms, the code I have put together below does not do partial pooling:

x2 = bf(X2 ~ 0 + X1) 
x3 = bf(X3 ~ 0 + X2)
y1 = bf(Y1 ~ 0 + X1)
y2 = bf(Y2 ~ 0 + Y1 + X2)
y3 = bf(Y3 ~ 0 + Y2 + X3)

fit1 = brm(data = XY,
              family = gaussian,
              formula = y1 + x2 + y2 + x3 + y3 + set_rescor(TRUE),
              prior = c(
                # x2 model
                prior(normal(0, 1), class = b, resp = X2),
                prior(exponential(1), class = sigma, resp = X2),
                
                # x3 model
                prior(normal(0, 1), class = b, resp = X3),
                prior(exponential(1), class = sigma, resp = X3),
                
                # y1 model
                prior(normal(0, 1), class = b, resp = Y1),
                prior(exponential(1), class = sigma, resp = Y1),
                
                # y2 model
                prior(normal(0, 1), class = b, resp = Y2),
                prior(exponential(1), class = sigma, resp = Y2),
                
                # y3 model
                prior(normal(0, 1), class = b, resp = Y3),
                prior(exponential(1), class = sigma, resp = Y3),
                
                # rho
                prior(lkj(2), class = rescor)
              )
)

Or, is there a way to manually create a grouping column with levels within the brmsformula? That way, I do not need to pass in multiple dataframes and a single dataframe will do. An example may be something like this (the non-working code below is just to demonstrate what I mean):

x2 = bf(X2 ~ 0 + (0 + X1 | level = "1", group = "alpha")) 
x3 = bf(X3 ~ 0 + (0 + X2 | level = "2", group = "alpha")) 
y1 = bf(Y1 ~ 0 + (0 + X1 | level = "1", group = "lambda"))
y2 = bf(Y2 ~ 0 + (0 + Y1 | level = "1", group = "beta") + (0 + X2 | level = "2", group = "lambda"))
y3 = bf(Y3 ~ 0 + (0 + Y2 | level = "2", group = "beta") + (0 + X3 | level = "3", group = "lambda"))

P.S. I am still new to brms and Bayesian modelling, so feel free to correct me or give feedback on my math notation and brms code.

Any help will be appreciated.

  • Operating System: 64-bit operating system, x64-based processor
  • brms Version: brms_2.15.0
2 Likes
#2

Hi,
sorry for not getting to you earlier, this is a good question!

Yes, there is! If I understand you correctly, then if we add a dummy column with a constant value then something like:

bf(X2 ~ 0 + (0 + X1 | gr(dummy, id  = "alpha")) +
bf(X3 ~ 0 + (0 + X2 | gr(dummy, id = "alpha")) + ...

would almost work. Specifying the same “id” value means that you build a shared correlation structure across all uses of the same ID. Problem is that dummy has only one level, but we will still estimate the sd for the dummy grouping, which is redundant. We can avoid this by setting the SD to a constant via the prior argument (this then serves the same role as putting a normal prior on the coefficient without grouping).

A quick simulated example:

dd <- data.frame(X1 = rnorm(10), X2 = rnorm(10), X3 = rnorm(10), dummy = "A")
brm(bf(X2 ~ 0 + ( 0 + X1 | gr(dummy, id = "alpha"))) + 
      bf(X3 ~ 0 + ( 0 + X1 | gr(dummy, id = "alpha"))) + 
      set_rescor(TRUE), data =dd)

And the output shows that we estimate the effects and their correlation + the residual sd and correlations but nothing more!

Group-Level Effects: 
~dummy (Number of levels: 1) 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(X2_X1)            3.00      0.00     3.00     3.00 1.00     4000     4000
sd(X3_X1)            3.00      0.00     3.00     3.00 1.00     4000     4000
cor(X2_X1,X3_X1)     0.00      0.67    -0.98     0.98 1.00     2044     1844

Family Specific Parameters: 
         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_X2     1.26      0.35     0.77     2.11 1.00     2462     1742
sigma_X3     1.15      0.31     0.71     1.92 1.00     2551     2305

Residual Correlations: 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(X2,X3)    -0.28      0.28    -0.76     0.33 1.00     2123     2074

Does that resolve your issue? Or am I missing what you need?

Best of luck with your model!

3 Likes
#3

Thank you so much Martin! Your reply has been super helpful!! :)

#4

@cao when you have a final model please post it here if possible. When people search they will fint it :)

1 Like