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:

- estimate the variance–covariance matrix between all variables by fitting a multivariate model with the MCMCglmm package, with all variables as response variables
- 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
- use the estimated matrices in the path analysis, that we perform using the sem package
- run each path models once for each of the 1000 estimated variance-covariance matrices
- extract the path coefficients and 95% confidence interval
- 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
- assess indirect effects by multiplying all coefficients on the focal path to obtain a path coefficient

Code:

```
#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,
nitt=1300000,
burnin=30000,
thin=100,
verbose = FALSE,
data = as.data.frame(data))
summary(mcmc1)
# 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)
```

- 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?