Specifying family in a PSEM run using brms

Interfaces brms
1

Hi! I am trying to run a piecewise structural equation model using five submodels and I was wondering how to get the complete model to run if the response variables in each submodels have different distributions (three are Gaussian, one is beta, and one is binomial). How can I specify family for each submodel in brm?

Any help and guidance is much appreciated! If there is a better way to code a PSEM in brms please let me know, I’m very new to the package!

Here is the code for each of my submodels and the complete model:

# Submodel 1 (gaussian)
hatch_mod <- bf(Hatch.Day ~ Supp.Factor + pc1.Fall +  pc1.Pre +  Breed.MT + (1|Nest.ID))

# Submodel 2 (gaussian)
brood_mod <- bf(Sibling ~ Supp.Factor + pc1.Fall +  pc1.Pre +  Breed.MT + Hatch.Day)

# Submodel 3 (beta)
diet_mod <- bf(Vert ~ Supp.Factor + pc1.Fall +  pc1.Pre +  Breed.MT + Hatch.Day + (1|Nest.ID))

# Submodel 4 (gaussian)
condition_mod <- bf(Condition ~ Vert + Supp.Factor + pc1.Fall +  pc1.Pre +  Breed.MT + Hatch.Day + Sibling + (1|Nest.ID))

# Submodel 5 (binomial)
survival_mod <- bf(Surv.Fall ~ Condition + Vert + Supp.Factor + pc1.Fall +  pc1.Pre +  Breed.MT + Hatch.Day + Sibling + (1|Nest.ID))

# Complete model
k_fit_brms <- brm(hatch_mod +
                    brood_mod + 
                    diet_mod +
                    condition_mod +
                    survival_mod +
                    set_rescor(FALSE), 
                  data=FinalDiet,
                  cores=2, chains = 2)
2

Are you using set_priors to set the priors for each model?

So my SEM looks like this

cw_mod <- bf(cw ~ fire_impact + flood_impact + cleared_impact + mean_depth_gw_neg +
               (1|site) + s(year, bs='gp', m=1) + (1|age))
will_mod <- bf(will ~ fire_impact + flood_impact + cleared_impact + mean_depth_gw_neg + 
                 (1|site) + s(year, bs='gp', m=1))
nmol_mod <- bf(nmol ~ fire_impact + flood_impact + cleared_impact + mean_depth_gw_neg + 
                 (1|site) + s(year, bs='gp', m=1))
ro_mod <- bf(ro ~ fire_impact + flood_impact + cleared_impact + mean_depth_gw_neg + 
               (1|site) + s(year, bs='gp', m=1))
sc_mod <- bf(sc ~ fire_impact + flood_impact + cleared_impact  + mean_depth_gw_neg + 
               (1|site) + s(year, bs='gp', m=1))
woody_mod <- bf(woody ~ fire_impact + flood_impact + cleared_impact + mean_depth_gw_neg + 
                  (1|site) + s(year, bs='gp', m=1))

k_fit_brms <- brm(cw_mod + will_mod + sc_mod + ro_mod + woody_mod + nmol_mod +
                    set_rescor(TRUE), 
                  data=clean_bemp_data, prior = prior2, family = "gaussian",
                  cores=6, chains = 4, iter = 3000,
                  control = list(adapt_delta=0.99, max_treedepth = 12))
3

With my priors set like this

prior2 <- c(set_prior("normal(-1, 1)", class = "b", coef= "fire_impact", resp="cw"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "flood_impact", resp="cw"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "cleared_impact", resp="cw"),
            set_prior("normal(0.5, 0.5)", class = "b", coef= "mean_depth_gw_neg",
                      resp="cw"),
            set_prior("normal(-1, 1)", class = "b", coef= "fire_impact", resp="will"),
            set_prior("normal(0.5, 0.5)", class = "b", coef= "flood_impact", resp="will"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "cleared_impact", resp="will"),
            set_prior("normal(0.5, 0.5)", class = "b", coef= "mean_depth_gw_neg",
                      resp="will"),
            set_prior("normal(-1, 1)", class = "b", coef= "fire_impact", resp="sc"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "flood_impact", resp="sc"),
            set_prior("normal(-1, 0.5)", class = "b", coef= "cleared_impact", resp="sc"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "mean_depth_gw_neg",
                      resp="sc"),
            set_prior("normal(-1, 1)", class = "b", coef= "fire_impact", resp="ro"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "flood_impact", resp="ro"),
            set_prior("normal(-1, 0.5)", class = "b", coef= "cleared_impact", resp="ro"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "mean_depth_gw_neg",
                      resp="ro"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "fire_impact", resp="nmol"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "flood_impact", resp="nmol"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "cleared_impact", resp="nmol"),
            set_prior("normal(0.5, 0.5)", class = "b", coef= "mean_depth_gw_neg",
                      resp="nmol"),
            set_prior("normal(1, 0.5)", class = "b", coef= "fire_impact", resp="woody"),
            set_prior("normal(-1, 0.5)", class = "b", coef= "flood_impact", resp="woody"),
            set_prior("student_t(7, 0, 2.5)", class = "b", coef= "cleared_impact", resp="woody"),
            set_prior("normal(-1, 0.5)", class = "b", coef= "mean_depth_gw_neg",
                      resp="woody"))
4

I all of my priors are uninformative.

Where I’m running into issues is that I’m not sure if I should be setting the family of brm() to gaussian when not all of the response variables in bf() are gaussian. One is beta and another is binomial.

5

Well the easy part first, you shouldn’t use uninformative piors. Spelling them out explicitly in the code is best practices.

Let me check to see what can be done with the response variables.

6

So in each bf() call you set the family there.

See here Brms: multivariate response, different distributions