Hello, I am using brms package to estimate multivariate models, each using different samples (i.e., the data observations used in the two are different).
I found that I could compare two or more coefficient from different models by using suest (seemingly unrelated estimation) and gsem (generalized structural equation modeling) functions in STATA.
I wonder if there is a similar options in brms.
Thanks.
Can you use an additional factor which indicates which sample the data come from and interact that with the other predictors in the model? Then everything’s in the same model and it should be easy to compare coefficients from the two data sets.
hi there. I don’t think that would work. It’s mainly because that different samples include same individuals; for example, my dv1 counts for free-service behavior while dv2 counts for paid-service behavior, among the same data panels.
Is there anybody who can help me?
Some more info about exactly what data you have would be helpful. But from what you’ve described so far it sounds like you could run a multivariate regression on your two DVs which includes random effects varying by participant. Within that model you will be able to compare the two sets of coefficients for the two DVs. You can correlate the random effects between the two DVs using the |ID| syntax as detailed here Estimating Multivariate Models with brms.
@andymilne Thank you so much for your kindness.
To be specific, my model is as follows :
model =
bf(Y1 ~ IV+(1|rand|individual)+factor(time), family=bernoulli(logit)) +
bf(Y2 ~ IV+(1|rand|individual)+factor(time), family=bernoulli(logit)) +
bf(DV1 | subset(Y1) ~ IV+(1|rand|individual)+factor(time), family=gaussian()) +
bf(DV2 | subset(Y1) ~ IV+(1|rand|individual)+factor(time), family=gaussian())
Prior_int = c(
set_prior(“normal(0,5)”, resp = "Y1, class = “Intercept”),
set_prior(“normal(0,5)”, resp = “Y2”, class = “Intercept”),
set_prior(“normal(0,5)”, resp = “DV1”, class = “Intercept”),
set_prior(“normal(0,5)”, resp = “DV2”, class = “Intercept”) )
brm (model + set_rescor(FALSE), data=data, […] , prior = c(Prior_int), inits=0)
, where Y1 and Y2 are dichotomous variables while DV1 and DV2 are continuous ones.
My concern is to compare the effect size of IV between Y1 and Y2, as well as D1 and D2.
I found ‘suest’ function in STATA do the similar things, but it does not exactly match with my object.
The formula looks reasonable (although I’ve not previously used the subset argument, so cannot comment on that).
This model presumably results in coefficients called something like Y1_IV and Y2_IV, and you want to compare those. Same with the coefficients DV1_IV and DV2_IV. The easiest way to compare coefficients is to use the hypothesis function:
hypothesis(model, c("Y1_IV - Y2_IV > 0", "DV1_IV - DV2_IV > 0"))
(or flipping the > to < depending on your hypothesized direction). This will give their estimated difference and CIs for that difference, along with an evidence ratio (posterior odds the difference is in the direction specified) and associated posterior probability. Hopefully I have understood what you are trying to do.
Thank you so much. Problem solved.
Hello again, I bring up a question related to the previous one that I have asked in this post. Considering my research objective, I need to compare the marginal effect of IV on Y1 and Y2, whether there is any difference in the impact size of the IV on two Ys which are non-linear. In this perspective, I wonder using the ‘hypothesis’ function can be an appropriate approach because as far as I know, it does not account for the differences in other results such as different intercept values.
I am not sure I fully understand your question. But note that, by default, hypothesis (for Y1 and Y2 in your model above) are done on the logit scale and all effects are linear on that scale.
If you want implied probabilities or differences of probabilities of “success” (or whatever 1 represents in Y1 and Y2), you can use the inv_logit_scaled transform within the hypothesis function and, within that, add the intercept and IV values you want. For example,
hypothesis(
model,
"(inv_logit_scaled(0.55 * Y1_IV) + 1.1) -
(inv_logit_scaled(0.33 * Y2_IV) + 2.2) > 0")
This would give the difference (and its CI) between the probability of Y1 “success” for IV = 0.55 with an intercept of 1.1, and the probability of Y2 “success” for IV = 0.33 with an intercept of 2.2. Obviously, you can plug in whichever IV and intercept values you need.
Edit – just to clarify that these tests do not take account of the random effects. If you wish to include those you are probably best off comparing predictions made with posterior_epred or posterior_predict.