Testing probability that a given situation occur

[EmiliaIOPAN]

Hi,

I am wondering, if it is possible to statistically compare the contribution of each source to the consumer diet from two different habitats. I run MixSIAR with few species (random factor) from vegetated and unvegetated habitat (fixed factor). I receive the estimated contribution of each sources for each species diet from two habitats but I would like to test the probability that the proportion of source X in the diet of consumer X from group X, is bigger than that in group Y. It was possible to test in SIAR, using this code:

Comparing the diets of the consumers in groups 4, 5 and 6 In this example i want to copmare the proportion of zostera in the diet across these 3 groups. These data are held in model1$output[,19] model1$output[,25] model1$output[,31] These data are depicted in the figure titled "Proportions by source: Zostera" To compare group 4 with group 5 i use the following code to deterimine which samples in the MCMC process were bigger in group 4 than group 5

test45 <- model1$output[,19] > model1$output[,25] the probability that the proportion of zostera in the diet of consumer group 4 is bigger than that in group 5 is then approximated by the proportion of samples that were bigger in group 4 than group 5 given by:

P45 <- sum(test45)/length(test45)

We can now do the exact same for groups 4 and 6 test46 <- model1$output[,19] > model1$output[,31] P46 <- sum(test46)/length(test46)

... or 5 & 6 test56 <- model1$output[,25] > model1$output[,31] P56 <- sum(test56)/length(test56)

... or 7 & 8 - which look very similar on the graph test78 <- model1$output[,37] > model1$output[,43] P78 <- sum(test78)/length(test78)

Thanks for any help.

Besr regards,

Emilia

[brianstock]

Hi Emilia,

Yes, you can do the same in MixSIAR, since it also generates posterior distributions for the contribution of each source (MCMC samples/chains). For examples on how to access the MCMC chains, see https://github.com/brianstock/MixSIAR/issues/30, https://github.com/brianstock/MixSIAR/issues/80, and https://github.com/brianstock/MixSIAR/issues/21. Also on pages 22-23 in the manual.

https://github.com/brianstock/MixSIAR/issues/21 shows how to calculate probabilities of interest from the MCMC chains (in wolves ex, probability that the Region 2 Deer diet proportion is greater than 0.7). In the same way as the SIAR example you show, you could calculate the probability that Deer contribute more to wolves diets in Region 1 than Region 2:

Run wolves example:

library(MixSIAR)
mixsiar.dir <- find.package("MixSIAR")
source(paste0(mixsiar.dir,"/example_scripts/mixsiar_script_wolves.R"))

Access posterior MCMC chains:

# p.fac1 indexed as [draws, region, source]
# median proportion Deer in Region 2 wolves diet
median(p.fac1[,2,1])
[1] 0.494242
# median proportion Deer in Region 1 wolves diet
median(p.fac1[,1,1])
[1] 0.7591761

This makes sense, we can see Region 1 wolves eat more Deer (manual p.20, also note I used "test" MCMC length here so the posteriors aren't pretty).

Now calculate probability that Deer contribute more to Region 1 diet than Region 2, as in SIAR ex above:

test12 <- p.fac1[,1,1] > p.fac1[,2,1]
P12 <- sum(test12)/length(test12)
P12
[1] 0.9193333