And now the full code…
# Calculate sumamry statistics for each group: TA, SEA and SEAc
group.ML <- groupMetricsML(siber.example)
print(group.ML)
1.1 1.2 1.3 1.4 2.1 2.2 2.3 2.4
TA 1.5555360 14.86029 84.20837 43.81102 4.552414 2.952808 4.434037 2.889293
SEA 0.5412062 5.35464 33.73813 15.87191 2.526144 1.666699 2.733655 1.873229
SEAc 0.5712732 5.65212 35.61247 16.75368 2.947168 1.944482 3.189264 2.185434So far these still just point-metrics that describe the width of the isotopic niche. That is, they are single numbers for each group, which means that we can’t compare one group to another in a statisical sense as we lack a measure of the uncertainty around each estimate. This is where we can use Bayesian Inference to quantify the error associated with fitting these ellipses to each group, that arises from both the number of samples we have, and also their distribution.
Essentially, what the MCMC algorithm does is generate a distribution of covariance matrices that to a greater or lesser extent (in terms of likelihood) describe the observed data. It does so, as is the general case in Bayesian Inference, by combing the prior probability with the likelihood of the data for a given covariance matrix.
SIBER uses the jags package to fit the Bayesian model and so we need to specify the parameters of the simulation run, including: run length, burn-in period, number of chains etc…
# options for running jags
parms <- list()
parms$n.iter <- 2 * 10^4 # number of iterations to run the model for
parms$n.burnin <- 1 * 10^3 # discard the first set of values
parms$n.thin <- 10 # thin the posterior by this many
parms$n.chains <- 2 # run this many chains
# define the priors
priors <- list()
priors$R <- 1 * diag(2)
priors$k <- 2
priors$tau.mu <- 1.0E-3
# fit the ellipses which uses an Inverse Wishart prior
# on the covariance matrix Sigma, and a vague normal prior on the
# means. Fitting is via the JAGS method.
ellipses.posterior <- siberMVN(siber.example, parms, priors)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 20
Unobserved stochastic nodes: 3
Total graph size: 35
Initializing model
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 20
Unobserved stochastic nodes: 3
Total graph size: 35
Initializing model
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 20
Unobserved stochastic nodes: 3
Total graph size: 35
Initializing model
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 20
Unobserved stochastic nodes: 3
Total graph size: 35
Initializing model
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 8
Unobserved stochastic nodes: 3
Total graph size: 23
Initializing model
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 8
Unobserved stochastic nodes: 3
Total graph size: 23
Initializing model
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 8
Unobserved stochastic nodes: 3
Total graph size: 23
Initializing model
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 8
Unobserved stochastic nodes: 3
Total graph size: 23
Initializing model
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%What we end up with is a range of ellipses that could explain the data, with more of them clustered around the most likely solution. However, one cannot simply take an average across these covariance matrices, as there are strict mathematical properties that must be maintained. The result of this is that it is not possible to plot a mean, median or modal Bayesian Standard Ellipse; instead we must calculate each one of the ellipse’s area, and then present summary statistics of this derived measurement. SIBER contains a function that will automatically loop over all the groups and do this.
The plots below represent the posterior distribution of the SEA_B fitted to each of the 4 groups in our dataset.
print(SEA.B.credibles)
$`1.1`
[,1] [,2]
99% 0.3468344 1.2578147
95% 0.4145694 1.0642662
50% 0.5792948 0.7926569
$`1.2`
[,1] [,2]
99% 2.695211 9.862806
95% 3.215251 8.356366
50% 4.456144 6.063702
$`1.3`
[,1] [,2]
99% 16.84132 62.43690
95% 20.42559 52.80049
50% 28.22531 38.37596
$`1.4`
[,1] [,2]
99% 7.540360 29.54346
95% 9.392405 24.40328
50% 13.330239 18.19873
$`2.1`
[,1] [,2]
99% 0.8228697 6.608253
95% 1.0868164 5.079182
50% 1.7760559 2.942953
$`2.2`
[,1] [,2]
99% 0.5279448 4.588260
95% 0.7206480 3.420861
50% 1.1969013 1.986865
$`2.3`
[,1] [,2]
99% 0.9058082 7.480594
95% 1.1680925 5.575338
50% 1.9357711 3.194633
$`2.4`
[,1] [,2]
99% 0.6003830 5.136584
95% 0.8128201 3.801781
50% 1.3142534 2.212923In order to test whether one group’s ellipse is smaller or larger than another, we can simply calculate the probability that its posterior distribution is smaller (or larger). This is acheived by comparing each pair of posterior draws for both groups, and dtermining which is smaller in magnitude. We then find the proportion of draws that are smaller, and this is a direct proxy for the probability that one group’s posterior distribution (of ellipse size in this case) is smaller than the other.
Here, we first calculate the proportion, and hence probability, of the SEA.B for group 1 being smaller than the SEA.B for group 2.
Pg1.lt.g2 <- sum( SEA.B[,1] < SEA.B[,2] ) / nrow(SEA.B)
print(Pg1.lt.g2)
[1] 1So, in this case, all of the estimates for groups 1’s ellipse are smaller than for group 2; although we could probably guess at this given that there appears to be no overlap between then 95% credible intervals of the two groups (see the figure above).
Then we can do exactly the same for groups 1 and 3.
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuUGcxLmx0LmczIDwtIHN1bSggU0VBLkJbLDFdIDwgU0VBLkJbLDNdICkgLyBucm93KFNFQS5CKVxucHJpbnQoUGcxLmx0LmczIClcbmBgYCJ9 -->
```r
Pg1.lt.g3 <- sum( SEA.B[,1] < SEA.B[,3] ) / nrow(SEA.B)
print(Pg1.lt.g3 )
```
<!-- rnb-source-end -->```r
Pg1.lt.g3 <- sum( SEA.B[,1] < SEA.B[,3] ) / nrow(SEA.B)
print(Pg1.lt.g3 )
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiWzFdIDFcbiJ9 -->
[1] 1
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
And then for the other pairings:
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVVR2MxTG14MExtYzNJRHd0SUhOMWJTZ2dVMFZCTGtKYkxEVmRJRHdnVTBWQkxrSmJMRGRkSUNrZ0x5QnVjbTkzS0ZORlFTNUNLVnh1Y0hKcGJuUW9VR2MxTG14MExtYzNLVnh1WEc1Z1lHQWlmUT09IC0tPlxuXG5gYGByXG5QZzUubHQuZzcgPC0gc3VtKCBTRUEuQlssNV0gPCBTRUEuQlssN10gKSAvIG5yb3coU0VBLkIpXG5wcmludChQZzUubHQuZzcpXG5cbmBgYFxuXG48IS0tIHJuYi1zb3VyY2UtZW5kIC0tPlxuIn0= -->
Pg5.lt.g7 <- sum( SEA.B[,5] < SEA.B[,7] ) / nrow(SEA.B)
print(Pg5.lt.g7)
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVZR0JnY2x4dVVHYzFMbXgwTG1jM0lEd3RJSE4xYlNnZ1UwVkJMa0piTERWZElEd2dVMFZCTGtKYkxEZGRJQ2tnTHlCdWNtOTNLRk5GUVM1Q0tWeHVjSEpwYm5Rb1VHYzFMbXgwTG1jM0tWeHVYRzVnWUdCY2JtQmdZQ0o5IC0tPlxuXG5gYGByXG5gYGByXG5QZzUubHQuZzcgPC0gc3VtKCBTRUEuQlssNV0gPCBTRUEuQlssN10gKSAvIG5yb3coU0VBLkIpXG5wcmludChQZzUubHQuZzcpXG5cbmBgYFxuYGBgXG5cbjwhLS0gcm5iLXNvdXJjZS1lbmQgLS0+XG4ifQ== -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuUGc1Lmx0Lmc3IDwtIHN1bSggU0VBLkJbLDVdIDwgU0VBLkJbLDddICkgLyBucm93KFNFQS5CKVxucHJpbnQoUGc1Lmx0Lmc3KVxuXG5gYGBcbmBgYCJ9 -->
```r
```r
Pg5.lt.g7 <- sum( SEA.B[,5] < SEA.B[,7] ) / nrow(SEA.B)
print(Pg5.lt.g7)
```
```
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiWzFdIDAuNTYzMjVcbiJ9 -->
```
[1] 0.56325
```
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
***
## Overlap Between Ellipses
One can calculate the overlap between two (or more) ellipses. In the first instance, this overlap is simply the area, in units of per mil squared, contained by the shape that lies within the overlapping region. This overlap is most easily calculated by using the SEAc of each ellipse.
The overlap between the SEAc for groups 3 and 4 in Community 1 is given by:
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVYRzV2ZG1WeWJHRndMa2N6TGtjMElEd3RJRzFoZUV4cGEwOTJaWEpzWVhBb1hDSXhMak5jSWl3Z1hDSXhMalJjSWl3Z2MybGlaWEl1WlhoaGJYQnNaU3dnY0NBOUlEQXVPVFVzSUc0Z1BURXdNQ2xjYmx4dVlHQmdJbjA9IC0tPlxuXG5gYGByXG5cbm92ZXJsYXAuRzMuRzQgPC0gbWF4TGlrT3ZlcmxhcChcXDEuM1xcLCBcXDEuNFxcLCBzaWJlci5leGFtcGxlLCBwID0gMC45NSwgbiA9MTAwKVxuXG5gYGBcblxuPCEtLSBybmItc291cmNlLWVuZCAtLT5cbiJ9 -->
overlap.G3.G4 <- maxLikOverlap(\1.3\, \1.4\, siber.example, p = 0.95, n =100)
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVZR0JnY2x4dVhHNXZkbVZ5YkdGd0xrY3pMa2MwSUR3dElHMWhlRXhwYTA5MlpYSnNZWEFvWEZ3eExqTmNYQ3dnWEZ3eExqUmNYQ3dnYzJsaVpYSXVaWGhoYlhCc1pTd2djQ0E5SURBdU9UVXNJRzRnUFRFd01DbGNibHh1WUdCZ1hHNWdZR0FpZlE9PSAtLT5cblxuYGBgclxuYGBgclxuXG5vdmVybGFwLkczLkc0IDwtIG1heExpa092ZXJsYXAoXFwxLjNcXCwgXFwxLjRcXCwgc2liZXIuZXhhbXBsZSwgcCA9IDAuOTUsIG4gPTEwMClcblxuYGBgXG5gYGBcblxuPCEtLSBybmItc291cmNlLWVuZCAtLT5cbiJ9 -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuXG5vdmVybGFwLkczLkc0IDwtIG1heExpa092ZXJsYXAoXFwxLjNcXCwgXFwxLjRcXCwgc2liZXIuZXhhbXBsZSwgcCA9IDAuOTUsIG4gPTEwMClcblxuYGBgXG5gYGAifQ== -->
```r
```r
overlap.G3.G4 <- maxLikOverlap(\1.3\, \1.4\, siber.example, p = 0.95, n =100)
```
```
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
And the overlap between SEAc of groups 1.2 and 2.1 is given by:
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHViM1psY214aGNDNURNVWN5TGtNeVJ6RWdQQzBnYldGNFRHbHJUM1psY214aGNDaGNJakV1TWx3aUxDQmNJakl1TVZ3aUxDQnphV0psY2k1bGVHRnRjR3hsTENCd0lEMGdNQzQ1TlN3Z2JpQTlJREV3TUNsY2JtQmdZQ0o5IC0tPlxuXG5gYGByXG5vdmVybGFwLkMxRzIuQzJHMSA8LSBtYXhMaWtPdmVybGFwKFxcMS4yXFwsIFxcMi4xXFwsIHNpYmVyLmV4YW1wbGUsIHAgPSAwLjk1LCBuID0gMTAwKVxuYGBgXG5cbjwhLS0gcm5iLXNvdXJjZS1lbmQgLS0+XG4ifQ== -->
overlap.C1G2.C2G1 <- maxLikOverlap(\1.2\, \2.1\, siber.example, p = 0.95, n = 100)
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVZR0JnY2x4dWIzWmxjbXhoY0M1RE1VY3lMa015UnpFZ1BDMGdiV0Y0VEdsclQzWmxjbXhoY0NoY1hERXVNbHhjTENCY1hESXVNVnhjTENCemFXSmxjaTVsZUdGdGNHeGxMQ0J3SUQwZ01DNDVOU3dnYmlBOUlERXdNQ2xjYm1CZ1lGeHVZR0JnSW4wPSAtLT5cblxuYGBgclxuYGBgclxub3ZlcmxhcC5DMUcyLkMyRzEgPC0gbWF4TGlrT3ZlcmxhcChcXDEuMlxcLCBcXDIuMVxcLCBzaWJlci5leGFtcGxlLCBwID0gMC45NSwgbiA9IDEwMClcbmBgYFxuYGBgXG5cbjwhLS0gcm5iLXNvdXJjZS1lbmQgLS0+XG4ifQ== -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxub3ZlcmxhcC5DMUcyLkMyRzEgPC0gbWF4TGlrT3ZlcmxhcChcXDEuMlxcLCBcXDIuMVxcLCBzaWJlci5leGFtcGxlLCBwID0gMC45NSwgbiA9IDEwMClcbmBgYFxuYGBgIn0= -->
```r
```r
overlap.C1G2.C2G1 <- maxLikOverlap(\1.2\, \2.1\, siber.example, p = 0.95, n = 100)
```
```
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
One might then wish to calculate the proportion overlap; athough one then runs into a choice as to what the demoninator will be in the equation. You could for instance calculate the proportion of A that overlaps with B, the proporiton of B that overlaps with A, or the proportion of A and B that overlap with each other.
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVjSEp2Y0M1dlppNW1hWEp6ZENBOExTQmhjeTV1ZFcxbGNtbGpLRzkyWlhKc1lYQXVSek11UnpSYlhDSnZkbVZ5YkdGd1hDSmRJQzhnYjNabGNteGhjQzVITXk1SE5GdGNJbUZ5WldFdU1Wd2lYU2xjYm5CeWFXNTBLSEJ5YjNBdWIyWXVabWx5YzNRcFhHNWdZR0FpZlE9PSAtLT5cblxuYGBgclxucHJvcC5vZi5maXJzdCA8LSBhcy5udW1lcmljKG92ZXJsYXAuRzMuRzRbXFxvdmVybGFwXFxdIC8gb3ZlcmxhcC5HMy5HNFtcXGFyZWEuMVxcXSlcbnByaW50KHByb3Aub2YuZmlyc3QpXG5gYGBcblxuPCEtLSBybmItc291cmNlLWVuZCAtLT5cbiJ9 -->
prop.of.first <- as.numeric(overlap.G3.G4[\overlap\] / overlap.G3.G4[\area.1\])
print(prop.of.first)
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVZR0JnY2x4dWNISnZjQzV2Wmk1bWFYSnpkQ0E4TFNCaGN5NXVkVzFsY21saktHOTJaWEpzWVhBdVJ6TXVSelJiWEZ4dmRtVnliR0Z3WEZ4ZElDOGdiM1psY214aGNDNUhNeTVITkZ0Y1hHRnlaV0V1TVZ4Y1hTbGNibkJ5YVc1MEtIQnliM0F1YjJZdVptbHljM1FwWEc1Z1lHQmNibUJnWUNKOSAtLT5cblxuYGBgclxuYGBgclxucHJvcC5vZi5maXJzdCA8LSBhcy5udW1lcmljKG92ZXJsYXAuRzMuRzRbXFxvdmVybGFwXFxdIC8gb3ZlcmxhcC5HMy5HNFtcXGFyZWEuMVxcXSlcbnByaW50KHByb3Aub2YuZmlyc3QpXG5gYGBcbmBgYFxuXG48IS0tIHJuYi1zb3VyY2UtZW5kIC0tPlxuIn0= -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucHJvcC5vZi5maXJzdCA8LSBhcy5udW1lcmljKG92ZXJsYXAuRzMuRzRbXFxvdmVybGFwXFxdIC8gb3ZlcmxhcC5HMy5HNFtcXGFyZWEuMVxcXSlcbnByaW50KHByb3Aub2YuZmlyc3QpXG5gYGBcbmBgYCJ9 -->
```r
```r
prop.of.first <- as.numeric(overlap.G3.G4[\overlap\] / overlap.G3.G4[\area.1\])
print(prop.of.first)
```
```
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiWzFdIDAuMTUwNjA0OVxuXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVjSEp2Y0M1dlppNXpaV052Ym1RZ1BDMGdZWE11Ym5WdFpYSnBZeWh2ZG1WeWJHRndMa2N6TGtjMFcxd2liM1psY214aGNGd2lYU0F2SUc5MlpYSnNZWEF1UnpNdVJ6UmJYQ0poY21WaExqSmNJbDBwWEc1d2NtbHVkQ2h3Y205d0xtOW1Mbk5sWTI5dVpDbGNibUJnWUNKOSAtLT5cblxuYGBgclxucHJvcC5vZi5zZWNvbmQgPC0gYXMubnVtZXJpYyhvdmVybGFwLkczLkc0W1xcb3ZlcmxhcFxcXSAvIG92ZXJsYXAuRzMuRzRbXFxhcmVhLjJcXF0pXG5wcmludChwcm9wLm9mLnNlY29uZClcbmBgYFxuXG48IS0tIHJuYi1zb3VyY2UtZW5kIC0tPlxuIn0= -->
[1] 0.1506049
prop.of.second <- as.numeric(overlap.G3.G4[\overlap\] / overlap.G3.G4[\area.2\])
print(prop.of.second)
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVZR0JnY2x4dWNISnZjQzV2Wmk1elpXTnZibVFnUEMwZ1lYTXViblZ0WlhKcFl5aHZkbVZ5YkdGd0xrY3pMa2MwVzF4Y2IzWmxjbXhoY0Z4Y1hTQXZJRzkyWlhKc1lYQXVSek11UnpSYlhGeGhjbVZoTGpKY1hGMHBYRzV3Y21sdWRDaHdjbTl3TG05bUxuTmxZMjl1WkNsY2JtQmdZRnh1WUdCZ0luMD0gLS0+XG5cbmBgYHJcbmBgYHJcbnByb3Aub2Yuc2Vjb25kIDwtIGFzLm51bWVyaWMob3ZlcmxhcC5HMy5HNFtcXG92ZXJsYXBcXF0gLyBvdmVybGFwLkczLkc0W1xcYXJlYS4yXFxdKVxucHJpbnQocHJvcC5vZi5zZWNvbmQpXG5gYGBcbmBgYFxuXG48IS0tIHJuYi1zb3VyY2UtZW5kIC0tPlxuIn0= -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucHJvcC5vZi5zZWNvbmQgPC0gYXMubnVtZXJpYyhvdmVybGFwLkczLkc0W1xcb3ZlcmxhcFxcXSAvIG92ZXJsYXAuRzMuRzRbXFxhcmVhLjJcXF0pXG5wcmludChwcm9wLm9mLnNlY29uZClcbmBgYFxuYGBgIn0= -->
```r
```r
prop.of.second <- as.numeric(overlap.G3.G4[\overlap\] / overlap.G3.G4[\area.2\])
print(prop.of.second)
```
```
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiWzFdIDAuMzIwMTMzM1xuXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVjSEp2Y0M1dlppNWliM1JvSUR3dElHRnpMbTUxYldWeWFXTW9iM1psY214aGNDNUhNeTVITkZ0Y0ltOTJaWEpzWVhCY0lsMGdMeUFvYjNabGNteGhjQzVITXk1SE5GdGNJbUZ5WldFdU1Wd2lYU0FySUc5MlpYSnNZWEF1UnpNdVJ6UmJYQ0poY21WaExqSmNJbDBwS1Z4dWNISnBiblFvY0hKdmNDNXZaaTVpYjNSb0tWeHVZR0JnSW4wPSAtLT5cblxuYGBgclxucHJvcC5vZi5ib3RoIDwtIGFzLm51bWVyaWMob3ZlcmxhcC5HMy5HNFtcXG92ZXJsYXBcXF0gLyAob3ZlcmxhcC5HMy5HNFtcXGFyZWEuMVxcXSArIG92ZXJsYXAuRzMuRzRbXFxhcmVhLjJcXF0pKVxucHJpbnQocHJvcC5vZi5ib3RoKVxuYGBgXG5cbjwhLS0gcm5iLXNvdXJjZS1lbmQgLS0+XG4ifQ== -->
[1] 0.3201333
prop.of.both <- as.numeric(overlap.G3.G4[\overlap\] / (overlap.G3.G4[\area.1\] + overlap.G3.G4[\area.2\]))
print(prop.of.both)
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVZR0JnY2x4dWNISnZjQzV2Wmk1aWIzUm9JRHd0SUdGekxtNTFiV1Z5YVdNb2IzWmxjbXhoY0M1SE15NUhORnRjWEc5MlpYSnNZWEJjWEYwZ0x5QW9iM1psY214aGNDNUhNeTVITkZ0Y1hHRnlaV0V1TVZ4Y1hTQXJJRzkyWlhKc1lYQXVSek11UnpSYlhGeGhjbVZoTGpKY1hGMHBLVnh1Y0hKcGJuUW9jSEp2Y0M1dlppNWliM1JvS1Z4dVlHQmdYRzVnWUdBaWZRPT0gLS0+XG5cbmBgYHJcbmBgYHJcbnByb3Aub2YuYm90aCA8LSBhcy5udW1lcmljKG92ZXJsYXAuRzMuRzRbXFxvdmVybGFwXFxdIC8gKG92ZXJsYXAuRzMuRzRbXFxhcmVhLjFcXF0gKyBvdmVybGFwLkczLkc0W1xcYXJlYS4yXFxdKSlcbnByaW50KHByb3Aub2YuYm90aClcbmBgYFxuYGBgXG5cbjwhLS0gcm5iLXNvdXJjZS1lbmQgLS0+XG4ifQ== -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucHJvcC5vZi5ib3RoIDwtIGFzLm51bWVyaWMob3ZlcmxhcC5HMy5HNFtcXG92ZXJsYXBcXF0gLyAob3ZlcmxhcC5HMy5HNFtcXGFyZWEuMVxcXSArIG92ZXJsYXAuRzMuRzRbXFxhcmVhLjJcXF0pKVxucHJpbnQocHJvcC5vZi5ib3RoKVxuYGBgXG5gYGAifQ== -->
```r
```r
prop.of.both <- as.numeric(overlap.G3.G4[\overlap\] / (overlap.G3.G4[\area.1\] + overlap.G3.G4[\area.2\]))
print(prop.of.both)
```
```
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiWzFdIDAuMTAyNDIxM1xuXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVjSEp2Y0M1dlppNWliM1JvTG14bGMzTXViM1psY214aGNDQThMU0JoY3k1dWRXMWxjbWxqS0c5MlpYSnNZWEF1UnpNdVJ6UmJYQ0p2ZG1WeWJHRndYQ0pkSUM4Z0tHOTJaWEpzWVhBdVJ6TXVSelJiWENKaGNtVmhMakZjSWwwZ0t5QnZkbVZ5YkdGd0xrY3pMa2MwVzF3aVlYSmxZUzR5WENKZElDMGdiM1psY214aGNDNUhNeTVITkZ0Y0ltOTJaWEpzWVhCY0lsMHBLVnh1Y0hKcGJuUW9jSEp2Y0M1dlppNWliM1JvTG14bGMzTXViM1psY214aGNDbGNibUJnWUNKOSAtLT5cblxuYGBgclxucHJvcC5vZi5ib3RoLmxlc3Mub3ZlcmxhcCA8LSBhcy5udW1lcmljKG92ZXJsYXAuRzMuRzRbXFxvdmVybGFwXFxdIC8gKG92ZXJsYXAuRzMuRzRbXFxhcmVhLjFcXF0gKyBvdmVybGFwLkczLkc0W1xcYXJlYS4yXFxdIC0gb3ZlcmxhcC5HMy5HNFtcXG92ZXJsYXBcXF0pKVxucHJpbnQocHJvcC5vZi5ib3RoLmxlc3Mub3ZlcmxhcClcbmBgYFxuXG48IS0tIHJuYi1zb3VyY2UtZW5kIC0tPlxuIn0= -->
[1] 0.1024213
prop.of.both.less.overlap <- as.numeric(overlap.G3.G4[\overlap\] / (overlap.G3.G4[\area.1\] + overlap.G3.G4[\area.2\] - overlap.G3.G4[\overlap\]))
print(prop.of.both.less.overlap)
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVZR0JnY2x4dWNISnZjQzV2Wmk1aWIzUm9MbXhsYzNNdWIzWmxjbXhoY0NBOExTQmhjeTV1ZFcxbGNtbGpLRzkyWlhKc1lYQXVSek11UnpSYlhGeHZkbVZ5YkdGd1hGeGRJQzhnS0c5MlpYSnNZWEF1UnpNdVJ6UmJYRnhoY21WaExqRmNYRjBnS3lCdmRtVnliR0Z3TGtjekxrYzBXMXhjWVhKbFlTNHlYRnhkSUMwZ2IzWmxjbXhoY0M1SE15NUhORnRjWEc5MlpYSnNZWEJjWEYwcEtWeHVjSEpwYm5Rb2NISnZjQzV2Wmk1aWIzUm9MbXhsYzNNdWIzWmxjbXhoY0NsY2JtQmdZRnh1WUdCZ0luMD0gLS0+XG5cbmBgYHJcbmBgYHJcbnByb3Aub2YuYm90aC5sZXNzLm92ZXJsYXAgPC0gYXMubnVtZXJpYyhvdmVybGFwLkczLkc0W1xcb3ZlcmxhcFxcXSAvIChvdmVybGFwLkczLkc0W1xcYXJlYS4xXFxdICsgb3ZlcmxhcC5HMy5HNFtcXGFyZWEuMlxcXSAtIG92ZXJsYXAuRzMuRzRbXFxvdmVybGFwXFxdKSlcbnByaW50KHByb3Aub2YuYm90aC5sZXNzLm92ZXJsYXApXG5gYGBcbmBgYFxuXG48IS0tIHJuYi1zb3VyY2UtZW5kIC0tPlxuIn0= -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucHJvcC5vZi5ib3RoLmxlc3Mub3ZlcmxhcCA8LSBhcy5udW1lcmljKG92ZXJsYXAuRzMuRzRbXFxvdmVybGFwXFxdIC8gKG92ZXJsYXAuRzMuRzRbXFxhcmVhLjFcXF0gKyBvdmVybGFwLkczLkc0W1xcYXJlYS4yXFxdIC0gb3ZlcmxhcC5HMy5HNFtcXG92ZXJsYXBcXF0pKVxucHJpbnQocHJvcC5vZi5ib3RoLmxlc3Mub3ZlcmxhcClcbmBgYFxuYGBgIn0= -->
```r
```r
prop.of.both.less.overlap <- as.numeric(overlap.G3.G4[\overlap\] / (overlap.G3.G4[\area.1\] + overlap.G3.G4[\area.2\] - overlap.G3.G4[\overlap\]))
print(prop.of.both.less.overlap)
```
```
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiWzFdIDAuMTE0MTA4NVxuIn0= -->
```
[1] 0.1141085
```
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
A problem with this simple overlap calculation is that it yields a point-estimate of overlap based on the maximum likelihood estimated SEA_c. One can instead calculate a distribution of overlap based on the posterior distirbutions of the fitted ellipses. It can be a bit slow to calculate this overlap, so you may want to drop the number of `draws` if your computer is slow.
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVZbUY1WlhNdWIzWmxjbXhoY0M1SE15NUhOQ0E4TFNCaVlYbGxjMmxoYms5MlpYSnNZWEFvWENJeExqTmNJaXdnWENJeExqUmNJaXhjYmlBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lHVnNiR2x3YzJWekxuQnZjM1JsY21sdmNpeGNiaUFnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUZ4dUlDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdaSEpoZDNNZ1BTQXhNREFzSUZ4dUlDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdjQzVwYm5SbGNuWmhiQ0E5SURBdU9UVXNYRzRnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNBZ0lDQWdJQ0FnSUNCdUlEMGdNell3S1Z4dWNISnBiblFvWW1GNVpYTXViM1psY214aGNDNUhNeTVITkNsY2JtQmdZQ0o5IC0tPlxuXG5gYGByXG5iYXllcy5vdmVybGFwLkczLkc0IDwtIGJheWVzaWFuT3ZlcmxhcChcIjEuM1wiLCBcIjEuNFwiLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZWxsaXBzZXMucG9zdGVyaW9yLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBkcmF3cyA9IDEwMCwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwLmludGVydmFsID0gMC45NSxcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG4gPSAzNjApXG5wcmludChiYXllcy5vdmVybGFwLkczLkc0KVxuYGBgXG5cbjwhLS0gcm5iLXNvdXJjZS1lbmQgLS0+XG4ifQ== -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYmF5ZXMub3ZlcmxhcC5HMy5HNCA8LSBiYXllc2lhbk92ZXJsYXAoXCIxLjNcIiwgXCIxLjRcIixcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGVsbGlwc2VzLnBvc3RlcmlvcixcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZHJhd3MgPSAxMDAsIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcC5pbnRlcnZhbCA9IDAuOTUsXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBuID0gMzYwKVxucHJpbnQoYmF5ZXMub3ZlcmxhcC5HMy5HNClcbmBgYCJ9 -->
```r
bayes.overlap.G3.G4 <- bayesianOverlap("1.3", "1.4",
ellipses.posterior,
draws = 100,
p.interval = 0.95,
n = 360)
print(bayes.overlap.G3.G4)
```
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiICAgICAgIGFyZWExICAgICBhcmVhMiAgIG92ZXJsYXBcbjEgICAyMjcuMTc1NCAgODcuODY3NTUgMjcuNzkzODczXG4yICAgMTczLjkwNDEgMTE3LjQ5NDUyIDI2LjA3MzE5OVxuMyAgIDIzMS44Mzc4ICA5My4wMzEyNCA0NS44NTE1MTVcbjQgICAxNzAuOTM0MiAgOTkuOTE3MDUgMTguNjM4MjU0XG41ICAgMTM0LjQwNTAgMTE0LjIxMDk1IDExLjU5MTMwMFxuNiAgIDE4NS4xMzMzICA3NC4xNDcxOSAyMi44MjI0MDZcbjcgICAxNDguNzE4NCAxMDEuNzA3NTMgMTYuNjc4NjM1XG44ICAgMTYxLjM0ODYgIDg3Ljg1Mjg1IDEyLjA5MzIxOVxuOSAgIDIyOC44MzMzICA2Ny40MDM2OCAyOS42OTgxNDNcbjEwICAxNjkuNzcxOCAxMTQuOTc5MDggMjQuODY0MDMzXG4xMSAgMTU1LjYzMjYgMTA1LjMwOTU0IDQyLjUxMjA0NVxuMTIgIDI1Ni43NjI1ICA4Ni44MTU1OSA0NC41OTA4MzVcbjEzICAxODcuODE1NCAxMDguNzgxNjIgMjYuNDU2NjMzXG4xNCAgMTUyLjA4MzEgMTAwLjY1MjU2IDE3Ljc3MTQ0OVxuMTUgIDIxNi41NDgyICA4NC4xNDc5MSAxNi40NDg4OTdcbjE2ICAyMzEuMDU0NiAxMjMuOTM5MDMgNTUuMDk0MTIzXG4xNyAgMjIyLjEwNDkgIDk5LjIzNzk2IDQ4LjU0NzQ1M1xuMTggIDMzOC43NTY0ICA5MS4yNTEzNyA1NS42OTY5ODBcbjE5ICAyMjAuMDUwOSAxMjUuNzY3ODYgMjIuNDg1NzM1XG4yMCAgMjAzLjc4ODIgIDg4LjgxNzM5IDIyLjc2Njg3M1xuMjEgIDIwOC44OTI5IDE5MS43NDE5NiA1Ni41NDUwNzZcbjIyICAyMjAuOTA5NSAgODcuNjgxNDIgNDUuMzcyOTc5XG4yMyAgMjAyLjE0MjMgMTA5LjYyMzUwIDMwLjA0MDAxOFxuMjQgIDE1Mi40NDk4IDExMS43OTcwOSAxOS42NTQ3ODVcbjI1ICAzMTMuODQ4NSAxMDAuODgwNTkgNjcuNDk4ODAwXG4yNiAgMTI1LjQzNzQgIDYyLjM5ODY4ICAxLjc4NzY4MlxuMjcgIDE5NS43NTkzICA5OS43MDUzNSAzMS40NjM4NDZcbjI4ICAyNjIuNDg4MiAgOTYuMDIzMTQgMzUuNTE0Nzk1XG4yOSAgNDMxLjYzMTIgIDkwLjUxNTY5IDQ5LjA3ODc5NFxuMzAgIDIzMy40MjM0IDEwMy4wNjA0MSA1MS43OTU1MjRcbjMxICAxNTEuNTEzOSAgOTguNDIwODQgIDUuNjgxMTgwXG4zMiAgMTY1LjAzNzYgIDkyLjM5MjY5IDE2LjUxNDI1M1xuMzMgIDE1OS41MDI4IDE0NS45MTE4NCAyOS4zMzAzNTNcbjM0ICAyNjMuNDIwNSAxMDUuMDkzODUgMzYuNTI4ODEwXG4zNSAgMTM5LjY4MjcgMTMyLjgxMzM0IDExLjA4NjAxNFxuMzYgIDIxNi4wMzk3IDEyMi43NjY5NCA0My4zODgzNzhcbjM3ICAxNzIuNjYwNCAgNzkuOTI2NjUgMTIuNTk3NjI3XG4zOCAgMjYxLjE1MjggMTQ2LjExNzYwIDUyLjM2NTg2M1xuMzkgIDE4NC44MjU1IDExMC4wOTI4NyAyNy43NTAyNzdcbjQwICAyMTYuMTIyMiAgODYuMjIyNDcgMjEuNDA0NzgxXG40MSAgMjk1LjkwNDQgIDkxLjQ0ODc3IDUyLjg2ODIyMlxuNDIgIDIyMS44NDg1IDE0My4xMDU5MyA0NC44NzA3MDhcbjQzICAyMTYuNTM1OSAgODEuMDE1MTggMjEuNDc0ODU0XG40NCAgMjc2LjMxMzUgIDg0Ljk1NDUxIDM2LjgwMTg5NlxuNDUgIDIxNC4wNTgyICA5NC4yMTkxNyAxOC4wODY4MjVcbjQ2ICAxMzIuOTAwMyAxMjAuNDc2NjUgMTcuMjg3NjQ4XG40NyAgMTcxLjc5MTUgIDcwLjc4NDMwIDEzLjA3MDM1N1xuNDggIDE2MS42NTE2IDEwMS4wNzEzOCAyMS42NDY0NDdcbjQ5ICAxNjAuMTg5OSAgNzUuOTU4OTggIDYuNDAyODg1XG41MCAgMjA0LjA1MTAgIDg0LjY1MTI0IDExLjE0NzU4MVxuNTEgIDE3MS4xMTI2ICA5NC43ODk3MyAzMi45MTA0OTRcbjUyICAyNDQuNDI1NyAgOTQuMDU3NjAgNDUuMDYxNjc3XG41MyAgMTcxLjgyMTUgIDk0LjU5MDc5IDIwLjk2MTAzNFxuNTQgIDE2OS4zOTIzIDEwMC44ODU0NSAyOC40NDc0NjlcbjU1ICAxOTYuNDIxOSAgNzIuMzYyMDEgMTcuMzYxNDAxXG41NiAgMTk0LjMzMzYgIDY5Ljg3MTM1IDI2LjU5NDE5OVxuNTcgIDM3Mi40NzIxIDExMy4yNzEzMiA2My42MDg3NzNcbjU4ICAyMzUuOTY1NCAxNDYuNjM4MzkgMzcuMzM3OTE4XG41OSAgMTkyLjQzNDggIDkzLjM5MzkyIDI0LjQ3ODg4OFxuNjAgIDE0MC42MDQ2IDEwMi45NTY1MSAgOC4zMDk1NDhcbjYxICAxMzYuMTE2MSAxMTIuODA1NjkgMjUuMTI3Mjc0XG42MiAgMjIwLjQ5MjMgIDg4LjkyNzgxIDE3LjgwODIwOFxuNjMgIDI1OC44MTUwIDE0MS4wMjEwNyA3Ni4yNTUyNjVcbjY0ICAyMDIuNzMxNiAgNjguMzMwMzQgMjMuNDM2MTAwXG42NSAgMjQzLjE5NjAgMTIzLjUzMjk1IDQzLjc4OTg5MlxuNjYgIDE5Ny45Nzk2IDE0Ny41NjA3OCA0NC4zNzI5MDVcbjY3ICAxNjEuNDE2MiAxMTIuMzI3MDQgMjUuOTM0MTI4XG42OCAgMTQxLjYzNzYgIDgxLjI4ODA5IDE5Ljc4ODAzMFxuNjkgIDE5NC41MjY3IDE1NS4zNjY1MyAzOS44MjQ5MzJcbjcwICAzMzYuMTU3NSAxMzYuMTIzMjEgNzYuNTY4OTA1XG43MSAgMTUwLjUxMTEgIDczLjcxOTE0ICA3LjUzODY4NlxuNzIgIDIzNy4yNjY0IDExNS4zMTA3MSA0NC45ODM1MjdcbjczICAxNzYuNDY5MCAgNTIuMDIzMzIgMTUuNTM2NjMyXG43NCAgMjA2LjExMDggMTAxLjQzNzI0IDE0LjI1OTExMlxuNzUgIDE5MC4xNDQ0IDExNC45NzE4NSAzMy42MzE5NTZcbjc2ICAyMjAuMzY3NCAgODEuNjk1NTMgMjcuOTIxMjg1XG43NyAgMTg5LjI3MDAgMTIxLjk2NzIzIDI3Ljc2NTAzOVxuNzggIDMxNC4xMTkzICA4OS4xMTQyNiA2MS4wODc3Mzhcbjc5ICAyMzcuNzY5MSAxMDcuMjAxMjQgMzkuMTI2MzA4XG44MCAgMjQyLjY5ODAgMTA3LjE1NDc0IDU1LjY0OTM0M1xuODEgIDIxNy44MzcxIDExNS43NDY5NyAzMS4xOTkzMTVcbjgyICAxODQuMjEzMyAxMTEuMjA2NjUgMTcuMDg4Nzg0XG44MyAgMjUxLjEzMzggMTA2LjYyODg3IDYyLjczNDk4N1xuODQgIDIxMS42NTgyIDExMy44MDQ4MCAzOC42MDQ3MjJcbjg1ICAyNTIuMzQ2OSAgNzkuNTAwMDYgNDUuNzE1MjY2XG44NiAgMjMwLjcwMTAgMTc2LjU3OTU1IDcwLjI4MjMzMlxuODcgIDI3Mi40MzI0ICA2My4yNDIxMyAyNi4yMzc5MzNcbjg4ICAyMDkuMTAwNSAxMDQuMjQ3MzIgNDMuODE1ODIzXG44OSAgMjI2LjMxOTYgMTA4LjM3MDc2IDMyLjk4MDkwNlxuOTAgIDE2OC4xNjc5IDEyMy41ODA3NSA0Mi40Mjc3MDZcbjkxICAyMjIuNTg5NCAgNjkuNTI2MjIgMzUuMTYxOTg3XG45MiAgMjMyLjQwMjIgIDc3LjE1MTgyIDI1LjY2Njg1MVxuOTMgIDE3Ny44NzQ1ICA5Ni43OTMyOSAxNC40MjU2NjNcbjk0ICAxNDYuNDQ5MyAxMDUuNzYxMjggMTEuODI2OTUxXG45NSAgMjMwLjMxNjAgIDg5LjMxNzYxIDEwLjQ1NDgxOVxuOTYgIDE4My4wMjM2ICA3NC43NTE0MiAyMC4zNDUwNTJcbjk3ICAxMjQuODMwMiAgOTYuNjA2NjQgMTEuMjAzNDExXG45OCAgMjU1LjQ4NDIgMTAwLjgyNzA4IDI0LjU2MjkyMFxuOTkgIDI5OC43ODU1ICA4NC4xMTA1NiAyOC44ODI5ODhcbjEwMCAxOTkuNTI2NyAgODYuNjY1NDkgMzEuMDQ2NTA4XG4ifQ== -->
```
area1 area2 overlap
1 227.1754 87.86755 27.793873
2 173.9041 117.49452 26.073199
3 231.8378 93.03124 45.851515
4 170.9342 99.91705 18.638254
5 134.4050 114.21095 11.591300
6 185.1333 74.14719 22.822406
7 148.7184 101.70753 16.678635
8 161.3486 87.85285 12.093219
9 228.8333 67.40368 29.698143
10 169.7718 114.97908 24.864033
11 155.6326 105.30954 42.512045
12 256.7625 86.81559 44.590835
13 187.8154 108.78162 26.456633
14 152.0831 100.65256 17.771449
15 216.5482 84.14791 16.448897
16 231.0546 123.93903 55.094123
17 222.1049 99.23796 48.547453
18 338.7564 91.25137 55.696980
19 220.0509 125.76786 22.485735
20 203.7882 88.81739 22.766873
21 208.8929 191.74196 56.545076
22 220.9095 87.68142 45.372979
23 202.1423 109.62350 30.040018
24 152.4498 111.79709 19.654785
25 313.8485 100.88059 67.498800
26 125.4374 62.39868 1.787682
27 195.7593 99.70535 31.463846
28 262.4882 96.02314 35.514795
29 431.6312 90.51569 49.078794
30 233.4234 103.06041 51.795524
31 151.5139 98.42084 5.681180
32 165.0376 92.39269 16.514253
33 159.5028 145.91184 29.330353
34 263.4205 105.09385 36.528810
35 139.6827 132.81334 11.086014
36 216.0397 122.76694 43.388378
37 172.6604 79.92665 12.597627
38 261.1528 146.11760 52.365863
39 184.8255 110.09287 27.750277
40 216.1222 86.22247 21.404781
41 295.9044 91.44877 52.868222
42 221.8485 143.10593 44.870708
43 216.5359 81.01518 21.474854
44 276.3135 84.95451 36.801896
45 214.0582 94.21917 18.086825
46 132.9003 120.47665 17.287648
47 171.7915 70.78430 13.070357
48 161.6516 101.07138 21.646447
49 160.1899 75.95898 6.402885
50 204.0510 84.65124 11.147581
51 171.1126 94.78973 32.910494
52 244.4257 94.05760 45.061677
53 171.8215 94.59079 20.961034
54 169.3923 100.88545 28.447469
55 196.4219 72.36201 17.361401
56 194.3336 69.87135 26.594199
57 372.4721 113.27132 63.608773
58 235.9654 146.63839 37.337918
59 192.4348 93.39392 24.478888
60 140.6046 102.95651 8.309548
61 136.1161 112.80569 25.127274
62 220.4923 88.92781 17.808208
63 258.8150 141.02107 76.255265
64 202.7316 68.33034 23.436100
65 243.1960 123.53295 43.789892
66 197.9796 147.56078 44.372905
67 161.4162 112.32704 25.934128
68 141.6376 81.28809 19.788030
69 194.5267 155.36653 39.824932
70 336.1575 136.12321 76.568905
71 150.5111 73.71914 7.538686
72 237.2664 115.31071 44.983527
73 176.4690 52.02332 15.536632
74 206.1108 101.43724 14.259112
75 190.1444 114.97185 33.631956
76 220.3674 81.69553 27.921285
77 189.2700 121.96723 27.765039
78 314.1193 89.11426 61.087738
79 237.7691 107.20124 39.126308
80 242.6980 107.15474 55.649343
81 217.8371 115.74697 31.199315
82 184.2133 111.20665 17.088784
83 251.1338 106.62887 62.734987
84 211.6582 113.80480 38.604722
85 252.3469 79.50006 45.715266
86 230.7010 176.57955 70.282332
87 272.4324 63.24213 26.237933
88 209.1005 104.24732 43.815823
89 226.3196 108.37076 32.980906
90 168.1679 123.58075 42.427706
91 222.5894 69.52622 35.161987
92 232.4022 77.15182 25.666851
93 177.8745 96.79329 14.425663
94 146.4493 105.76128 11.826951
95 230.3160 89.31761 10.454819
96 183.0236 74.75142 20.345052
97 124.8302 96.60664 11.203411
98 255.4842 100.82708 24.562920
99 298.7855 84.11056 28.882988
100 199.5267 86.66549 31.046508
```
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
And summarise the credible intervals of the Bayesian overlap output. Note that this code does not work well on the small number of posterior draws we are using for this basic example - for one it returns negative values which is not possible, but is arising as the smoother has not got enough information to stay close to or within the positive number range.
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVJeUJoYm1RZ2QyVWdZMkZ1SUdOaGJHTjFiR0YwWlNCMGFHVWdZMjl5Y21WemNHOXVaR2x1WnlCamNtVmthV0pzWlNCcGJuUmxjblpoYkhNZ2RYTnBibWRjYmlNZ2IzVnlJR052WkdVZ1puSnZiU0JoWW05MlpTQmhaMkZwYmx4dUl5QmpZV3hzSUhSdklHaGtjbU5rWlRwb1pISWdkWE5wYm1jZ2JHRndjR3g1S0NsY2JtOTJaWEpzWVhBdVkzSmxaR2xpYkdWeklEd3RJR3hoY0hCc2VTaGNiaUFnWVhNdVpHRjBZUzVtY21GdFpTaGlZWGxsY3k1dmRtVnliR0Z3TGtjekxrYzBLU3dnWEc0Z0lHWjFibU4wYVc5dUtIZ3NMaTR1S1h0MGJYQThMV2hrY21Oa1pUbzZhR1J5S0hncEpHaGtjbjBzWEc0Z0lIQnliMklnUFNCamNpNXdLVnh1WEc1d2NtbHVkQ2h2ZG1WeWJHRndMbU55WldScFlteGxjeWxjYm1CZ1lDSjkgLS0+XG5cbmBgYHJcbiMgYW5kIHdlIGNhbiBjYWxjdWxhdGUgdGhlIGNvcnJlc3BvbmRpbmcgY3JlZGlibGUgaW50ZXJ2YWxzIHVzaW5nXG4jIG91ciBjb2RlIGZyb20gYWJvdmUgYWdhaW5cbiMgY2FsbCB0byBoZHJjZGU6aGRyIHVzaW5nIGxhcHBseSgpXG5vdmVybGFwLmNyZWRpYmxlcyA8LSBsYXBwbHkoXG4gIGFzLmRhdGEuZnJhbWUoYmF5ZXMub3ZlcmxhcC5HMy5HNCksIFxuICBmdW5jdGlvbih4LC4uLil7dG1wPC1oZHJjZGU6Omhkcih4KSRoZHJ9LFxuICBwcm9iID0gY3IucClcblxucHJpbnQob3ZlcmxhcC5jcmVkaWJsZXMpXG5gYGBcblxuPCEtLSBybmItc291cmNlLWVuZCAtLT5cbiJ9 -->
# and we can calculate the corresponding credible intervals using
# our code from above again
# call to hdrcde:hdr using lapply()
overlap.credibles <- lapply(
as.data.frame(bayes.overlap.G3.G4),
function(x,...){tmp<-hdrcde::hdr(x)$hdr},
prob = cr.p)
print(overlap.credibles)````
```r
# and we can calculate the corresponding credible intervals using
# our code from above again
# call to hdrcde:hdr using lapply()
overlap.credibles <- lapply(
as.data.frame(bayes.overlap.G3.G4),
function(x,...){tmp<-hdrcde::hdr(x)$hdr},
prob = cr.p)
print(overlap.credibles)
<!-- rnb-source-end -->
<!-- rnb-output-end -->
<!-- rnb-output-begin eyJkYXRhIjoiJGFyZWExXG4gICAgICAgICBbLDFdICAgICBbLDJdXG45OSUgIDg2LjQ4ODY5IDMyNi43NzUyXG45NSUgMTE5LjQxNzM2IDI4Ny44MTA4XG41MCUgMTY3LjU4NDc2IDIzNC4wMzIzXG5cbiRhcmVhMlxuICAgICAgICBbLDFdICAgICBbLDJdXG45OSUgNDQuNzI0NDggMTU3LjQ0NDBcbjk1JSA1Ny43MDMxNiAxMzguNDEwMFxuNTAlIDgwLjU1ODU5IDEwOC4wODE4XG5cbiRvdmVybGFwXG4gICAgICAgICBbLDFdICAgICBbLDJdXG45OSUgLTkuMjY0NTI5IDgxLjE3OTQ5XG45NSUgIDEuMDE4Mjc3IDYwLjUxMjUwXG41MCUgMTcuMjM1NDY3IDM0LjIxMjExXG4ifQ== -->
$area1 [,1] [,2] 99% 86.48869 326.7752 95% 119.41736 287.8108 50% 167.58476 234.0323
$area2 [,1] [,2] 99% 44.72448 157.4440 95% 57.70316 138.4100 50% 80.55859 108.0818
$overlap [,1] [,2] 99% -9.264529 81.17949 95% 1.018277 60.51250 50% 17.235467 34.21211 ```