How to run a path analysis in a bayesian framework?

Hello all!

I am trying to run a path analysis to understand how reproductive output (i.e., number of offspring that survive until adulthood) is directly versus indirectly related to exploration and parental investment in a frog species.
I want to use a methodology that has been used in a few papers so far (e.g., Mutzel et al., 2013; Zhao et al., 2016), but I am stuck in some steps.
The method is the following:
  1. estimate the variance–covariance matrix between all variables by fitting a multivariate model with the MCMCglmm package, with all variables as response variables
  2. run the same model with four different priors: (1) inverse Wishart (V = diag(n), ν = n), (2) inverse gamma (V = diag(n), ν = b.002), (3) flat prior (V = diag(n), ν = 1.002) and (4) parameter expanded (V = diag(n), ν = c), with n = number of response variables in the model, b = n ¬- 1 and c = n + 1
  3. use the estimated matrices in the path analysis, that we perform using the sem package
  4. run each path models once for each of the 1000 estimated variance-covariance matrices
  5. extract the path coefficients and 95% confidence interval
  6. when the 95% confidence intervals only slightly overlapped zero, check how many times the estimate was positive or negative, to get an equivalent of a p-value
  7. assess indirect effects by multiplying all coefficients on the focal path to obtain a path coefficient
Here is the code I wrote so far:
#In a frog species, I expect the reproductive output to be directly influenced by the number of times
an individual deposit tadpoles in a water body (nb.deposition), the number of different water bodies
an individual uses to deposit its tadpoles (nb.pools), the rapidity of the individual to discover new
water bodies (week.first.discovery). I also expect these 3 variables to be directly influenced by
the individual’s level of exploration. Therefore, exploration should indirectly influence the reproductive output.

#preparation of the four priors

prior1 = list(R = list(V = diag(5), nu = 5))

prior2 = list(R = list(V = diag(5), nu = 5.002))

prior3 = list(R = list(V = diag(5)*10-6, nu = 6))

prior4 = list(R = list(V = diag(5), nu = 6))

#estimation of the variance-covariance matrix between the five variables (I will repeat the same model with the four different priors)

mcmc1 <- MCMCglmm(cbind(output, nb.deposition, nb.pools, week.first.discovery, exploration) ~ trait-1,
                          rcov =~ us(trait-1):units,
                          family = c("poisson","poisson","poisson","poisson","gaussian"),
                          prior = prior1,
                          verbose = FALSE,
                          data =

# path analysis

sem.mod <- readLines(con = textConnection(
  "output ->  nb.deposition, lam1,  NA
output ->   nb.pools, lam2,  NA
output ->   week.first.discovery, lam3, NA
nb.deposition ->   exploration, lam4,  NA
nb.pools ->   exploration, lam5,  NA
week.first.discovery ->   exploration, lam6,  NA
nb.deposition <-> nb.deposition, e1,    NA
nb.pools <-> nb.pools, e2,    NA
week.first.discovery <-> week.first.discovery, e3,    NA
exploration <-> exploration, e4,    NA
output <-> output, NA,    1
nb.deposition <-> nb.deposition, NA,    1
nb.pools <-> nb.pools, NA,    1
week.first.discovery <-> week.first.discovery, NA,    1
cat(sem.mod, file = 'sem.mod.spec', sep = '\n')
sem.mod1<- specifyModel(file = 'sem.mod.spec')
sem.ind <- sem(sem.mod1, “variance-covariance matrices”, N)
So my questions are:
  • How can I get the variance-covariance matrices estimated in the MCMC and run each path models once for each of the 1000 estimated variance-covariance matrices?
  • How do I extract the path coefficients and 95% confidence interval?