Niche Overlap Analysis Using nicheROVER and Stable Isotope Data

Author

Affiliation

Department of Biology, University of Hamburg

Published

September 18, 2024

Abstract

Stable isotope analysis (SIA) provides numerous information while assessing species’ ecological niche, the R package nicheROVER offers an insightful analysis using stable isotope data.

Part of this note is taken from the course Community Ecology offered by University of Hamburg for M.Sc. Biology program and led by Dr. Paula Cabral Eterovick.

The cover image is credited to Museums Victoria on Unsplash.

R Setup

# setup
library(tidyverse)
library(nicheROVER)

Niche Overlap

In Ecological Niche, we talked about the basic ideas of ecological niches and niche overlap. Understanding the niche overlap patterns, both intraspecific and interspecific, can help us better understand many topics, such as:

Evolutionary:

Sexual Dimorphism by studying differences in resource utilization between males and females of the same species.
Speciation Process by studying niche partitioning between sister taxa or subspecies.

Conservation:

Impact of Invasive Species by examining niche overlap between native and invasive species.
Ecological Guild by understanding niche overlap within a group sharing the same or similar functional traits (e.g. pollinator communities).

However, comparing niches among multiple species is challenging, as I mentioned earlier, since ecological data are almost always multidimensional.

`nicheROVER`

The package nicheROVER provides a new probabilistic method for quantifying n-dimensional ecological niche and niche overlap (Swanson et al. 2015). The script demonstrated here is written by Martin Lysy, Heidi Swanson (see here), and is also used in their publication.

fish is a fish stable isotope dataset included in the package, containing values for three stable isotopes (N-15, C-13, S-34) measured in the muscle tissue of four Arctic fish species (ARCS, BDWF, LKWF, LSCS). More details can be found from the documentation (?fish) or in the publication.

Table 1: fish dataset

ABCDEFGHIJ0123456789

species <fct>	D15N <dbl>	D13C <dbl>	D34S <dbl>
ARCS	12.18	-23.89	13.07
ARCS	11.19	-23.79	13.77
ARCS	12.05	-22.88	14.03
ARCS	11.25	-24.98	14.18
ARCS	12.08	-25.82	11.84
ARCS	11.54	-25.69	11.23
ARCS	12.64	-22.72	12.46
ARCS	11.25	-24.79	14.62
ARCS	10.94	-26.99	11.26
ARCS	13.18	-22.47	11.74

What is SIA?

Stable isotope analysis (SIA) is a powerful tool in ecological niche analysis. Many physicochemical and biochemical processes are sensitive to differences in the dissociation energies of molecules, which often depend on the mass of the elements from which these molecules are made. Thus, the isotopic composition of many materials (expressed as $δ$ -values, indicating the amount ratio of certain isotope to a standard value), including the tissues of organisms, often contains a label of the process that created it (Newsome et al. 2007). These “signatures” could provide information about their ecological roles and interactions.

Table 2: Common isotope systems and expected patterns in

δ

-values, subsetted from (Newsome et al. 2007)

Gradient	Isotope system	High $δ$ -values	Low $δ$ -values
Trophic level	$δ^{13} C$ / $δ^{15} N$	High levels	Low levels
C₃-C₄ Vegetation	$δ^{13} C$	C₄ plants	C₃ plants
Marine-terrestrial	$δ^{15} N$ / $δ^{13} C$ / $δ^{34} S$	High levels	Low levels

For more information about isotopic and atomic ecology

Newsome et al. (2007)
Michel and Lepoint (2021)
Moyo et al. (2021)

Bayesian Workflow

Bayesian theorem

Van De Schoot et al. (2021) described the foundational idea of the Bayesian statistics as follows:

Bayesian statistics is an approach to data analysis based on Bayes’ theorem, where available knowledge about parameters in a statistical model is updated with the information in observed data. The background knowledge is expressed as a prior distribution and combined with observational data in the form of a likelihood function to determine the posterior distribution.

This means we make an initial estimation of our parameters (prior distribution), then observe data and update the model (posterior distribution). Bayes’ theorem, which this approach is based on, is mathematically stated as following:

$P (A | B) = \frac{P (B | A) P (A)}{P (B)}$

where,

$A$ and $B$ are events and $P (B) \neq 0$ ,
$P (A | B)$ : the probability of event $A$ occurring given that $B$ is true. It is also called the posterior probability (briefly posterior) of $A$ given $B$ .
$P (B | A)$ is also a conditional probability: the probability of event $B$ occurring given that $A$ is true. It can also be interpreted as the likelihood of $A$ given a fixed $B$ because $P (B | A) = L (A | B)$ . ( $L$ : likelihood function)
$P (A)$ and $P (B)$ are the probabilities of observing $A$ and $B$ respectively without any given conditions; they are known as the prior probability and marginal probability.

Why Bayesian Statistics?

Bayesian statistics is a different framework and approach on statistical analysis. Both classical and Bayesian inference, while trying to model data, have own advantages and disadvantages. Bayesian inference is generally employed whenever uncertainty comes into play:

There is also a practical advantage of the Bayesian approach, which is that all its inferences are probabilistic and thus can be represented by random simulations. For this reason, whenever we want to summarize uncertainty in estimation beyond simple confidence intervals, and whenever we want to use regression models for predictions, we go Bayesian.

― Gelman, Hill, and Vehtari (2024), P. 16

Swanson et al. (2015) pointed out a viewpoint from the traditional niche overlap analysis that might introduce problems: the geometric approach doesn’t account for uncertainty and also lacks direction. This is well illustrated by Fig. 1 in the publication:

Figure 1: Comparison between geometric and probabilistic isotope overlap metrics. Isotope x is distributed in species A and B (x-axis); dashed lines show means. The 95% niche region (solid lines) of B is completely contained within that of A, such that the geometric overlap of B onto A is 100%. However, the probability that a randomly chosen member from A falls into the isotopic niche of B is only 35%. (Swanson et al. 2015)

By inviting the variance of the niche probability distribution, it now also differentiates the overlap in directions, which means the probability of a species A falling into species B’s niche is not equal to that of species B falling into species A’s niche. Since this perspective starts to consider probability in the models, a Bayesian workflow will provide more reliable and robust analysis.

Random Vector and Multivariate Normal Distribution

Given an individual statistical unit, in our case, the niche described by three isotope ratios, a random vector is a group of variables that describes this mathematical system. Which is again, in our case, the three values of the isotope ratio. The distribution of each component of the random vector together forms a multivariate distribution.

When dealing with univariate (one-dimensional) analysis, we often assume that the variable is normally distributed. We do likewise even with multivariate analysis. The multivariate normal distribution is a generalization of the univariate normal distribution into higher dimensions, but it’s not as simple as all components of the random vector being normal. It also requires every linear combination of arbitrary two components in the system to be normally distributed (bivariate normal distribution), which is not always true even if every single variable is normal on its own. This is where the covariance matrix comes into play. Normal distribution is described by two variables: mean and standard deviation (sd), with the latter being simply the square root of the variance.

Covariance and covariance matrix

Covariance is a measure of the joint variability of two variables. The sign of the covariance shows the tendency, i.e. one of the three possible linear relationships between two variables, which are:

negative linear relationship, covariance is < 0
no linear relationship, covariance is around 0
positive linear relationship, covariance is > 0

The variance is a special case of covariance in which the two variables are identical. Covariance itself is sensitive to the scale of the data, and it does not provide information of on the extent to which the two variables correlate with each other and how well the linear correlation fits the data. All these characteristics make covariance difficult to interpret, but it has become an important stepping stone for many subsequent analysis, such as PCA, correlation analysis, etc.

A covariance matrix is a square matrix giving the covariance between each pair of elements, which is always symmetric, and the diagonal elements represent the variance of each variable. The covariance matrix describes the joint behavior of every linear combination of two variables, analogue to univariate normal distribution, a covariance matrix ( $Σ$ ) and a mean ( $μ$ ) are used to described a multivariate normal distribution.

For more information about covariance

StatQuest on YouTube (simply and clearly explained)

Conjugate Prior

Jackson et al. (2011) suggested that the normal distribution is generally appropriate when dealing with stable isotope ratios, and this gives us a way to simplify the analysis by taking advantage of the conjugate prior.

A conjungate prior is a certain prior distribution that, when applied to a certain likelihood function, results in a posterior distribution from the same distribution family. For example, in univaraite analysis, a normal prior and a normal likelihood function results a normal posterior. In our case, the conjugate prior for the multivariate normal is the normal-inverse-wishart (NIW) distribution. It can be considered in two parts: the mean follows a normal distribution and the covariance follows an inverse-Wishart distribution. NIW is commonly used as the prior distribution when modeling a multivariate normal distribution where the means and the covariances are unknown.

The probability of a random observation of a isotope ratio ( $x$ ) belongs to species $s$ is following a normal distribution.

$P (x | species = s) ∽ N (μ_{s}, Σ_{s})$

Monte Carlo with `niw.post()`

All the prior knowledge we just went through, from basic Bayesian inference to the NIW distribution, is to understand the very two lines of code that implemented the Bayesian workflow into the analysis.

nsamples <- 1000
fish.par <- tapply(1:nrow(fish), fish$species,
            function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4]))

niw.post() uses a data matrix (X) as observation data, calculate the posterior distribution, and then draws a set of means and covariance matrices of a possible multivariate distribution, which includes the prior and our observation data. The main idea is that because analytically working with posterior multivariate distribution is too impractical, so we employ the Monte Carlo approach to have a view of the posterior. We use tapply() to separately apply 1000 times of niw.post() on each fish species.

This results in a list of 4 lists (species) of 2 matrices (mu and Sigma), where the mu matrix has a dimension of 1000 (sampling times) x 3 (number of variables), and the Sigma matrix has a dimension of 1000 (sampling times) x 3 x 3 (one slice = one covariance matrix).

This particular data structure is critical for later application of niche.par.plot(), see ?nichr.par.plot for more.

See the result

str(fish.par)

List of 4
 $ ARCS:List of 2
  ..$ mu   : num [1:1000, 1:3] 12.6 12.7 12.6 12.5 12.6 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  ..$ Sigma: num [1:3, 1:3, 1:1000] 0.601 0.369 0.345 0.369 1.146 ...
  .. ..- attr(*, "dimnames")=List of 3
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  .. .. ..$ : NULL
 $ BDWF:List of 2
  ..$ mu   : num [1:1000, 1:3] 9.15 9.28 9.48 9.41 9.28 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  ..$ Sigma: num [1:3, 1:3, 1:1000] 1.251 -0.313 -0.415 -0.313 4.629 ...
  .. ..- attr(*, "dimnames")=List of 3
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  .. .. ..$ : NULL
 $ LKWF:List of 2
  ..$ mu   : num [1:1000, 1:3] 11 11 11.1 11.1 11.1 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  ..$ Sigma: num [1:3, 1:3, 1:1000] 0.916 0.106 0.73 0.106 2.305 ...
  .. ..- attr(*, "dimnames")=List of 3
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  .. .. ..$ : NULL
 $ LSCS:List of 2
  ..$ mu   : num [1:1000, 1:3] 11.8 11.6 11.6 11.6 11.9 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  ..$ Sigma: num [1:3, 1:3, 1:1000] 0.581 -0.0583 0.1468 -0.0583 1.6501 ...
  .. ..- attr(*, "dimnames")=List of 3
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  .. .. ..$ : chr [1:3] "D15N" "D13C" "D34S"
  .. .. ..$ : NULL
 - attr(*, "dim")= int 4
 - attr(*, "dimnames")=List of 1
  ..$ : chr [1:4] "ARCS" "BDWF" "LKWF" "LSCS"

Niche Parameters

niche.par.plot() allows us to have a overview on Monte Carlo simulation of the posterior multivariate distribution of our $μ$ and $Σ$ . It shows the probability density curve of each variable. Use plot.index, plot.mu or plot.Sigma to select which $μ$ or $Σ$ to plot, or simply use plot.mu = TRUE, plot.Sigma = TRUE to display all. See details in the document.

clrs <- c('#348E8C', '#732C02', '#BF0404', '#F29F05') # colors for each species

# mu1 (del15N), mu2 (del13C), and Sigma12 
par(mar = c(4, 5, .5, .1)+.1, oma = c(0,0,0,0), pty="s", mfrow = c(1,3))
niche.par.plot(fish.par, col = clrs, plot.index = 1)
niche.par.plot(fish.par, col = clrs, plot.index = 2)
niche.par.plot(fish.par, col = clrs, plot.index = 1:2)
legend("topleft", legend = names(fish.par), fill = clrs)

Figure 2: The distribution of $μ_{1}$ ( $δ^{15} N)$ and $μ_{2}$ ( $δ^{13} C)$ , as well as the $Σ_{12}$ is displayed here.

2-D Ellipse Projection

To deal with n-dimensional niche regions (N_R), Swanson et al. (2015) reduced them into 2-dimensional space using probabilistic projection, which is to differentiate from geometric projection. A 2-dimensional probabilistic projection of an n-dimensional N_R is simply the N_R calculated by the projected 2 axes (Jackson et al. 2011; Swanson et al. 2015). Fig. 2 from Swanson et al. (2015) clearly demonstrates and explains the different outcomes.

Figure 3: A three-dimensional 95% niche region for normally distributed isotopes. The two ellipses are projections of the niche region (NR) onto isotope axes U and V. The black, hatched ellipse is the geometric projection of NR onto this plane, which is the shadow obtained by shining a light onto the three-dimensional ellipsoid along the W-axis (in blue). This is the largest of the two-dimensional ellipse “slices” taken parallel to the UV-plane. The red ellipse is the 95% niche region constructed directly from isotopes U and V. Thus, the geometric projection of the ellipsoid covers less than 95% of randomly selected individuals in two-dimensional isotope space. (Swanson et al. 2015)

For the sake of clarity, we only use 10 simulations here. niche.plot() generates three types of plot:

Univariate distribution density curves, i.e. smoothed histogram.
Bivariate scatter plots.
2-D projection of N_R using simulation data, displayed in ellipse plots.

From the first two plot types, we could confirm the assumption that the isotope ratios are generally normal. The ellipse plots give us a better view of each species’ N_R and the overlap situation.

nsamples <- 10
fish.par <- tapply(1:nrow(fish), fish$species,
            function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4]))

# format data for plotting function
fish.data <- tapply(1:nrow(fish), fish$species, function(ii) X = fish[ii,2:4])

niche.plot(niche.par = fish.par, niche.data = fish.data, pfrac = .05,
           iso.names = expression(delta^{15}*N, delta^{13}*C, delta^{34}*S),
           col = clrs, xlab = expression("Isotope Ratio (per mil)"))

It’s a bit overwhelming plot. Take time to understand it.

Estimate Niche Overlap

With overlap() we can have a more accurate and robust inference of overlapping among the species. The accuracy of the estimation can be improved by increasing nprob. The generated plot shows pairwise niche overlaps (now we have directions!), which displays the probability distribution of species A (along rows) fallen inside the niche of species B (along columns), with indication of the mean and 95% credible interval.

Be careful not to confuse two “95%” here - one means that we use 95% of the niche region to calculate overlap, which is specified in overlap() with alpha = 0.95; the other is the 95% credible interval of the probability distribution displayed in each plot. There is a difference to spot between plots using 95% or 99% of the niche region.

nsamples <- 1000
fish.par <- tapply(1:nrow(fish), fish$species,
            function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4]))

# 95% niche region
over.stat95 <- overlap(fish.par, nreps = nsamples, nprob = 1000, alpha = 0.95)
overlap.plot(over.stat95, col = clrs, mean.cred.col = "#024873", equal.axis = TRUE, xlab = "Overlap Probability (%) -- Niche Region Size: 95%")

Figure 5: Niche overlap using 95% of the niche region.

# 99% niche region
over.stat99 <- overlap(fish.par, nreps = nsamples, nprob = 1000, alpha = 0.99)
overlap.plot(over.stat99, col = clrs, mean.cred.col = "#024873", equal.axis = TRUE, xlab = "Overlap Probability (%) -- Niche Region Size: 99%")

Figure 6: Niche overlap using 99% of the niche region.

By printing out the returned value from overlap(), numeric results can also be shown.

See detailed results

We could calculate the estimation for both alpha = 0.99 and 0.95 at the same time.

over.stat.2lvl <- overlap(fish.par, nreps = nsamples, nprob = 1000, alpha = c(0.95, 0.99))
over.mean.2lvl <- apply(over.stat.2lvl, c(1:2, 4), mean)*100

# Here are the mean of each overlap probability distribution
round(over.mean.2lvl, 2)

, , alpha = 95%

         Species B
Species A  ARCS  BDWF  LKWF  LSCS
     ARCS    NA 10.64 65.25 80.95
     BDWF  0.31    NA 25.72  4.50
     LKWF  7.35 78.18    NA 52.68
     LSCS 37.77 51.57 88.46    NA

, , alpha = 99%

         Species B
Species A  ARCS  BDWF  LKWF  LSCS
     ARCS    NA 32.33 86.64 91.71
     BDWF  0.79    NA 41.23  8.93
     LKWF 11.64 92.30    NA 70.00
     LSCS 50.48 80.34 97.02    NA

# Here are the upper and lower bounds of the 95% credible interval
over.cred <- apply(over.stat.2lvl*100, c(1:2, 4), quantile, prob = c(.025, .975), na.rm = TRUE)
round(over.cred)

, , Species B = ARCS, alpha = 95%

       Species A
        ARCS BDWF LKWF LSCS
  2.5%    NA    0    3   27
  97.5%   NA    1   12   49

, , Species B = BDWF, alpha = 95%

       Species A
        ARCS BDWF LKWF LSCS
  2.5%     1   NA   61   25
  97.5%   29   NA   92   79

, , Species B = LKWF, alpha = 95%

       Species A
        ARCS BDWF LKWF LSCS
  2.5%    38   15   NA   75
  97.5%   88   37   NA   97

, , Species B = LSCS, alpha = 95%

       Species A
        ARCS BDWF LKWF LSCS
  2.5%    66    2   38   NA
  97.5%   93    8   69   NA

, , Species B = ARCS, alpha = 99%

       Species A
        ARCS BDWF LKWF LSCS
  2.5%    NA    0    6   39
  97.5%   NA    2   18   64

, , Species B = BDWF, alpha = 99%

       Species A
        ARCS BDWF LKWF LSCS
  2.5%     9   NA   80   55
  97.5%   65   NA   99   97

, , Species B = LKWF, alpha = 99%

       Species A
        ARCS BDWF LKWF LSCS
  2.5%    66   27   NA   91
  97.5%   98   55   NA  100

, , Species B = LSCS, alpha = 99%

       Species A
        ARCS BDWF LKWF LSCS
  2.5%    82    4   55   NA
  97.5%   98   15   85   NA

Penguins

Let’s run the niche overlap analysis on another dataset. There’s a datasets of 3 Pygoscelis species nesting along the Palmer Archipelago near Palmer Station, containing structural measurements, two isotope ratios, and basic metadata. Isotope data was collected from blood samples, and additional records of whether the sampled individuals have completed their clutch were documented.

adeliae <- readRDS("biology-research/community-ecology/assets/data/219.5.rds") %>% mutate(Species = "P. adeliae")
papua <- readRDS("biology-research/community-ecology/assets/data/220.7.rds") %>% mutate(Species = "P. papua")
antarcticus <- readRDS("biology-research/community-ecology/assets/data/221.8.rds") %>% mutate(Species = "P. antarcticus")

penguin <- bind_rows(adeliae, papua, antarcticus)

comp.penguin <- penguin %>% select(c(1:9, 13:14, 17, 10:12, 15:16)) %>% filter(complete.cases(.))

paged_table(comp.penguin)

Table 3: penguin dataset

ABCDEFGHIJ0123456789

studyName <fct>	Sample.Number <int>	Species <chr>	Region <fct>	Island <fct>	Stage <fct>	Individual.ID <fct>	Clutch.Completion <fct>	Date.Egg <date>	Body.Mass <int>
PAL0708	2	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N1A2	Yes	2007-11-11	3800
PAL0708	3	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N2A1	Yes	2007-11-16	3250
PAL0708	5	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N3A1	Yes	2007-11-16	3450
PAL0708	6	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N3A2	Yes	2007-11-16	3650
PAL0708	7	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N4A1	No	2007-11-15	3625
PAL0708	8	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N4A2	No	2007-11-15	4675
PAL0708	10	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N5A2	Yes	2007-11-09	4250
PAL0708	11	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N6A1	Yes	2007-11-09	3300
PAL0708	15	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N8A1	Yes	2007-11-16	4400
PAL0708	17	P. adeliae	Anvers	Torgersen	Adult, 1 Egg Stage	N9A1	Yes	2007-11-12	3450

To consider potential sexual dimorphisms within species, we will also differentiate the observations by sex.

comp.penguin <- comp.penguin %>% 
    filter(Sex != "" & Sex != ".") %>%
    mutate(Species = factor(Species), Sex = factor(case_when(
        Sex == "FEMALE" ~ "F",
        Sex == "MALE" ~ "M"
    ))) %>% 
    mutate(cat = factor(paste(Species, Sex, sep = " - ")))

# Take a look at the n in each group
summary(comp.penguin$cat)

    P. adeliae - F     P. adeliae - M P. antarcticus - F P. antarcticus - M 
                71                 68                 34                 33 
      P. papua - F       P. papua - M 
                58                 60

As mentioned in Table 2, $δ^{13} C$ and $δ^{15} N$ can be used to assess the trophic level of the species. We first analyze the niche using only these two values:

Code

# Run the analysis only with isotope data
nsamples <- 10
pen.par <- tapply(
    1:nrow(comp.penguin), 
    comp.penguin$cat, 
    function(ii) niw.post(nsamples = nsamples, X = comp.penguin[ii, 16:17]))
    
# format data for plotting function
pen.data <- tapply(1:nrow(comp.penguin), comp.penguin$cat, function(ii) X = comp.penguin[ii, 16:17])

clrs <- c("#023059", "#0477BF", "#8C3E7F", "#8872A6", "#F24130", "#F2A71B")

par(mfcol = (c(6, 6)))
niche.plot(niche.par = pen.par, niche.data = pen.data, pfrac = .05,
           iso.names = c(expression(delta^{15}*N), expression(delta^{13}*C)),
           col = clrs, xlab = expression("Isotope Ratio (per mil)"))

nsamples <- 1000
pen.par <- tapply(
    1:nrow(comp.penguin), 
    comp.penguin$cat, 
    function(ii) niw.post(nsamples = nsamples, X = comp.penguin[ii, 16:17]))

# 95% niche region
over.stat95 <- overlap(pen.par, nreps = nsamples, nprob = 1000, alpha = 0.95)
overlap.plot(over.stat95, col = clrs, mean.cred.col = "#024873", equal.axis = TRUE, xlab = "Overlap Probability (%) -- Niche Region Size: 95%", species.names = c("P. adeliae - F", "P.adeliae - M", "P. antarcticus - F", "P. antarcticus - M", "P. papua - F", "P. papua - M"))

Figure 7: No sexual difference in $δ^{13} C$ and $δ^{15} N$ signature is shown. *P. adeliae* has the biggest niche region.

Figure 8: Niche overlap probability distribution among species and sexes, only isotope data used.

In Figure 7 we saw no remarkable differences between sexes within each species. We also noticed that P. adeliae has the largest niche region and it overlaps with both the others. This is also shown in Figure 8.

Then, we combine isotope data with two structural measurements, culmen length and culmen depth, whose functionality should be strongly associated with their diet, to run the analysis again.

Code

# Run the analysis with culmen size data and isotope data
nsamples <- 10
pen.par <- tapply(
    1:nrow(comp.penguin), 
    comp.penguin$cat, 
    function(ii) niw.post(nsamples = nsamples, X = comp.penguin[ii, c(13:14, 16:17)]))
    
# format data for plotting function
pen.data <- tapply(1:nrow(comp.penguin), comp.penguin$cat, function(ii) X = comp.penguin[ii, c(13:14, 16:17)])

clrs <- c("#023059", "#0477BF", "#8C3E7F", "#8872A6", "#F24130", "#F2A71B")

par(mfcol = (c(6, 6)))
niche.plot(niche.par = pen.par, niche.data = pen.data, pfrac = .05,
           iso.names = c("Culmen length", "Culmen depth", expression(delta^{15}*N), expression(delta^{13}*C)),
           col = clrs, xlab = expression("Isotope Ratio (per mil)"))

nsamples <- 1000
pen.par <- tapply(
    1:nrow(comp.penguin), 
    comp.penguin$cat, 
    function(ii) niw.post(nsamples = nsamples, X = comp.penguin[ii, c(13:14, 16:17)]))

# 95% niche region
over.stat95 <- overlap(pen.par, nreps = nsamples, nprob = 1000, alpha = 0.95)
overlap.plot(over.stat95, col = clrs, mean.cred.col = "#024873", equal.axis = TRUE, xlab = "Overlap Probability (%) -- Niche Region Size: 95%", species.names = c("P. adeliae - F", "P.adeliae - M", "P. antarcticus - F", "P. antarcticus - M", "P. papua - F", "P. papua - M"))

Figure 9: Niche plot using culmen size and isotope data.

Figure 10: Niche overlap probability using culmen size data and isotope data.

In Figure 9 we notice that both in culmen length and depth, a notable sexual dimorphism is present, and this has formed different niche regions among the groups. The overlap in niche region defined by this data subset (Figure 10) shows, however, different results than those which only utilized isotope signatures (Figure 8). The overlap in all groups and directions between P. papua and P. adeliae is gone, suggesting that their diets have similar $δ^{13} C$ and $δ^{15} N$ signatures, but the prey is probably not the same, since they evolved different functional traits in bill structure. This is more obvious in Figure 9, there is no overlap in niche region described by culmen depth and length between the two species.

On the other hand, female individuals of P. antarcticus remain the only group which has a notable overlap with other groups, namely both sexes of P. adeliae. Gorman, Williams, and Fraser (2014), who conducted the study utilizing this data, suggested, regarding their $δ^{15} N$ signature, that the female P. antarcticus might rely more on a certain type of prey, which might then overlap with P. adeliea’s menu. This might also be a reason for the decreased breeding success.

Table 4: Percentage of breed pair which completed the whole clutch (i.e. 2 eggs for Pygoscelis species)

clutch <- comp.penguin %>% group_by(Species) %>% summarise(complete_clutch = round(sum(Clutch.Completion == "Yes")/sum(!is.na(Clutch.Completion))*100, 2))
clutch

# A tibble: 3 × 2
  Species        complete_clutch
  <fct>                    <dbl>
1 P. adeliae                90.6
2 P. antarcticus            79.1
3 P. papua                  94.1

And at last, these are our adorable guests:

P. adeliae by Dylan Shaw on Unsplash — *P. adeliae* by Dylan Shaw on Unsplash

P. papua by Australian Antarctic Program — *P. papua* by Australian Antarctic Program

By the way, if anyone ever wondered – yes, this is the exact same dataset they used to produce the palmerpenguins package. I found the data at EDI Data Portal and didn’t notice until I attached the images, because these three guys showing up at the same time is weirdly familiar to me. The dataset is therefore available by simply attaching the package. If you are interested in what more detailed results people have found with these data, you should check the paper from Gorman, Williams, and Fraser (2014).

Surprise!

Reference

Gelman, Andrew, Jennifer Hill, and Aki Vehtari. 2024. “Regression and Other Stories (Corrections up to 26 Oct 2023) Please Do Not Reproduce in Any Form Without Permission,” July. https://avehtari.github.io/ROS-Examples/.

Gorman, Kristen B, Tony D Williams, and William R Fraser. 2014. “Ecological Sexual Dimorphism and Environmental Variability Within a Community of Antarctic Penguins (Genus Pygoscelis).” PLOS ONE 9 (3). https://doi.org/10.1371/journal.pone.0090081.

Jackson, Andrew L., Richard Inger, Andrew C. Parnell, and Stuart Bearhop. 2011. “Comparing Isotopic Niche Widths Among and Within Communities: SIBER - Stable Isotope Bayesian Ellipses in R: Bayesian Isotopic Niche Metrics.” Journal of Animal Ecology 80 (3): 595–602. https://doi.org/10.1111/j.1365-2656.2011.01806.x.

LTER, Palmer Station Antarctica, and Kristen Gorman. 2020a. “Structural Size Measurements and Isotopic Signatures of Foraging Among Adult Male and Female Adélie Penguins (Pygoscelis Adeliae) Nesting Along the Palmer Archipelago Near Palmer Station, 2007-2009.” Environmental Data Initiative. https://doi.org/10.6073/PASTA/98B16D7D563F265CB52372C8CA99E60F.

———. 2020b. “Structural Size Measurements and Isotopic Signatures of Foraging Among Adult Male and Female Chinstrap Penguins (Pygoscelis Antarcticus) Nesting Along the Palmer Archipelago Near Palmer Station, 2007-2009.” Environmental Data Initiative. https://doi.org/10.6073/PASTA/CE9B4713BB8C065A8FCFD7F50BF30DDE.

———. 2020c. “Structural Size Measurements and Isotopic Signatures of Foraging Among Adult Male and Female Gentoo Penguins (Pygoscelis Papua) Nesting Along the Palmer Archipelago Near Palmer Station, 2007-2009.” Environmental Data Initiative. https://doi.org/10.6073/PASTA/9FC8F9B5A2FA28BDCA96516649B6599B.

Michel, Loïc N., and Gilles Lepoint. 2021. “Stable Isotopes as Descriptors of Ecological Niches.” Zenodo. https://doi.org/10.5281/ZENODO.5765530.

Moyo, Sydney, Hayat Bennadji, Danielle Laguaite, Anna A. Pérez-Umphrey, Allison M. Snider, Andrea Bonisoli-Alquati, Jill A. Olin, et al. 2021. “Stable Isotope Analyses Identify Trophic Niche Partitioning Between Sympatric Terrestrial Vertebrates in Coastal Saltmarshes with Differing Oiling Histories.” PeerJ 9 (July): e11392. https://doi.org/10.7717/peerj.11392.

Newsome, Seth D., Carlos Martinez Del Rio, Stuart Bearhop, and Donald L. Phillips. 2007. “A Niche for Isotopic Ecology.” Frontiers in Ecology and the Environment 5 (8): 429–36. https://doi.org/10.1890/060150.1.

Swanson, Heidi K., Martin Lysy, Michael Power, Ashley D. Stasko, Jim D. Johnson, and James D. Reist. 2015. “A New Probabilistic Method for Quantifying n -Dimensional Ecological Niches and Niche Overlap.” Ecology 96 (2): 318–24. https://doi.org/10.1890/14-0235.1.

Van De Schoot, Rens, Sarah Depaoli, Ruth King, Bianca Kramer, Kaspar Märtens, Mahlet G. Tadesse, Marina Vannucci, et al. 2021. “Bayesian Statistics and Modelling.” Nat Rev Methods Primers 1 (1): 1. https://doi.org/10.1038/s43586-020-00001-2.

R session and R packages

R session

R version 4.4.1 (2024-06-14) Platform: aarch64-apple-darwin20 Running under: macOS 15.2

Matrix products: default BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0

locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Berlin tzcode source: internal

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] palmerpenguins_0.1.1 nicheROVER_1.1.2 lubridate_1.9.3
[4] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[7] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
[10] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
[13] knitr_1.48 rmarkdown_2.28 pacman_0.5.1

Packages

pacman

Version: 0.5.1

Rinker TW, Kurkiewicz D (2018). pacman: Package Management for R. version 0.5.0, http://github.com/trinker/pacman.

rmarkdown

Version: 2.28

Allaire J, Xie Y, Dervieux C, McPherson J, Luraschi J, Ushey K, Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2024). rmarkdown: Dynamic Documents for R. R package version 2.28, https://github.com/rstudio/rmarkdown.

Xie Y, Allaire J, Grolemund G (2018). R Markdown: The Definitive Guide. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9781138359338, https://bookdown.org/yihui/rmarkdown.

Xie Y, Dervieux C, Riederer E (2020). R Markdown Cookbook. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9780367563837, https://bookdown.org/yihui/rmarkdown-cookbook.

knitr

Version: 1.48

Xie Y (2024). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.48, https://yihui.org/knitr/.

Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/.

Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.

tidyverse

Version: 2.0.0

Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.

nicheROVER

Version: 1.1.2

Lysy M, Stasko A, Swanson H (2023). nicheROVER: Niche Region and Niche Overlap Metrics for Multidimensional Ecological Niches. R package version 1.1.2, https://CRAN.R-project.org/package=nicheROVER.

palmerpenguins

Version: 0.1.1

Horst AM, Hill AP, Gorman KB (2020). palmerpenguins: Palmer Archipelago (Antarctica) penguin data. doi:10.5281/zenodo.3960218 https://doi.org/10.5281/zenodo.3960218, R package version 0.1.0, https://allisonhorst.github.io/palmerpenguins/.

Citation

BibTeX citation:

@online{hu2024,
  author = {Hu, Zhehao},
  title = {Niche {Overlap} {Analysis} {Using} {nicheROVER} and {Stable}
    {Isotope} {Data}},
  date = {2024-09-18},
  url = {https://zzzhehao.github.io/biology-research/community-ecology/nicheROVER-Niche-Overlap.html},
  langid = {en},
  abstract = {Stable isotope analysis (SIA) provides numerous
    information while assessing species’ ecological niche, the R package
    *nicheROVER* offers an insightful analysis using stable isotope
    data.}
}

For attribution, please cite this work as:

Hu, Zhehao. 2024. “Niche Overlap Analysis Using nicheROVER and Stable Isotope Data.” September 18, 2024. https://zzzhehao.github.io/biology-research/community-ecology/nicheROVER-Niche-Overlap.html.