---
title: "Usage of the niche metric functions (former 'resniche' package)"
author: 
  - "Miquel De Cáceres, CREAF"
  - "Oriol Lapiedra, CREAF"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: TRUE
bibliography: resniche.bib  
vignette: >
  %\VignetteIndexEntry{Usage of the niche metric functions (former 'resniche' package)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignettePackage{indicspecies}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Introduction

In this document we show how to use the functions described in @DeCaceres2011 by following an example of dietary preferences in pigeons belonging to two populations (Moià and Barcelona) in Catalonia (north east of Spain). We begin by loading the library and the data:
```{r}
library(indicspecies)
data(pigeons)
```
The data consists of three items: the resource use matrix for each pigeon population and the matrix of dissimilarities between pairs of resources (seeds). 
```{r}
ls()
```

## Resemblance between diet resources
Before starting any resource niche calculation, we can inspect the matrix of dissimilarities between the $r=6$ resources, $\mathbf{D}$:
```{r}
dfood
```
These can be summarized using a dendrogram. For example:

```{r, fig=TRUE, fig.height=5, fig.width=5}
plot(hclust(dfood, method="average"), h=-1, xlab="", 
     ylab="Distance", main="", sub="", ylim=c(0.2,1))

```

Some seeds are quite similar (e.g. popcorn, wheat or oats), whereas sunflower differs substantially from all other resources.

##  Resource niche analysis at the population level 

### Resource use of populations
We begin by showing the resource use data for pigeons of Barcelona and Moia, expressed as proportions (i.e., the vector $\mathbf{f}$ for each population):
```{r}
diet.pop.barcelona <- colSums(diet.barcelona)
round(diet.pop.barcelona/sum(diet.pop.barcelona), dig=3)
diet.pop.moia <- colSums(diet.moia)
round(diet.pop.moia/sum(diet.pop.moia), dig=3)
```

Whereas pigeons in Moià feed almost exclusively on wheat, those of Barcelona combine wheat and sunflower seeds.

### Niche breadth in populations
We will determine the resource niche breadth of each of the two populations as @DeCaceres2011: 

$B_{pop} = (1/2)\sum_{j=1}^r\sum_{k=1}^r{f_jf_kd_{jk}}$

We first conduct our calculations without taking into account the resemblance between resources (which is equal to stating that $d_{jk}=1$ in all cases):
```{r}
nichevar(P=diet.barcelona, mode="single")
nichevar(P=diet.moia, mode="single")
```

In general, we can say that the niche breadth of the population in Moià is smaller than the niche breadth of the population in Barcelona. If we repeat the same calculations with the matrix of resource resemblance, we realize that the niche breadth of both populations becomes smaller:
```{r}
popvar.barcelona <- nichevar(P=diet.barcelona, D=dfood, 
                            mode="single")
popvar.barcelona
popvar.moia <- nichevar(P=diet.moia, D=dfood, mode="single")
popvar.moia
```

The reason is that the first analysis was assuming that all resources were equally (and maximally) distinct, while the second analysis accounts for the similarity between some resources. Moreover, note that the niche breadth of Moià has decreased more than the niche breadth of Barcelona. This reflects that the resources being used by pigeons of Moià are more similar than the resources used by pigeons of Barcelona.


### Overlap between populations
We can now calculate the niche overlap between the two pigeon populations:
\[
O_{12} = \frac{\sum_{j=1}^r\sum_{k=1}^r{f_{1j}f_{2k}(1-d_{jk}^2)}}{\sqrt{\sum_{j=1}^r\sum_{k=1}^r{f_{1j}f_{1k}(1-d_{jk}^2)}\sum_{j=1}^r\sum_{k=1}^r{f_{2j}f_{2k}(1-d_{jk}^2)}}}
\]
Using function `nicheoverlap()`:
```{r}
nicheoverlap(P1=diet.barcelona, P2=diet.moia, mode="single")
nicheoverlap(P1=diet.barcelona, P2=diet.moia, mode="single", D = dfood)
```

If we include the resemblance between resources, the degree of overlap increases, for the same reason that we obtained smaller niche breadth statistics when resemblances were included.

## Resource niche analysis at the individual level
In this section, we perform a resource niche analysis at individual level. In particular, we are interested in assessing how much the resource niche of individuals differs from that of their corresponding population. For this, we need to calculate a measure of the degree of individual specialization. 

### Resource use of individuals
We begin by showing the resource use data for all `r nrow(diet.barcelona)` pigeons in the sample from the Barcelona population, expressed as proportions (i.e., matrix $\mathbf{F}$ for Barcelona):

```{r}
round(sweep(diet.barcelona, 1, FUN="/", 
            rowSums(diet.barcelona)), dig=3)
```

We see that most individuals in Barcelona feed on either sunflower or wheat, but there are some individuals (like pigeon 9) which prefer oats. Now we display the resource use data for the `r nrow(diet.moia)` pigeons representing the population in Moià (i.e., matrix $\mathbf{F}$ for Moià):
```{r}
round(sweep(diet.moia, 1, FUN="/", 
            rowSums(diet.moia)), dig=3)
```

Many pigeons from Moià feed on wheat seeds almost exclusively, but some of them have broader preferences. 

### Measuring the degree of individual specialisation
We begin our resource niche analysis by calculating the niche breadth of each individual in the population, $B_i$:
\[
B_i = (1/2)\sum_{j=1}^r\sum_{k=1}^r{f_{ij}f_{ik}d_{jk}}
\]
We calculate the values for the individuals of both populations:
```{r}
indvar.barcelona <- nichevar(P=diet.barcelona, D=dfood)
summary(indvar.barcelona)
indvar.moia <- nichevar(P=diet.moia, D=dfood)
summary(indvar.moia)
```

Most individuals have niche breadths that are smaller than the niche breadth of their corresponding population, although a few individuals in Moià have niche breadths larger than the population value. A niche breadth equal to zero indicates that only one resource is exploited. 


We can compare the niche breadth values of the two populations using a non-parametric test.
```{r}
wilcox.test(indvar.barcelona$B, indvar.moia$B)
```

The Wilcoxon test confirms that the niche breadth of pigeons in Barcelona is generally higher than that of pigeons in Moià, as we saw at the population level.

In order to calculate the degree of individual specialization, @Bolnick2002 defined WIC/TNW, i.e. the ratio between the within individual component (i.e. average niche width) and the total niche width of the population. Similarly we define the following specialization measure, that takes into account the resemblance between resources:
\[
S_{pop} = \frac{\sum_{i=n}^n{B_i}/n}{B_{pop}}
\]
where $B_i$ is the niche breadth of each individual, and $B_{pop}$ is the niche breadth of the population. Note that it is possible that $B_i$ values can be larger than $B_{pop}$. However, we do not expect the average of $B_i$ values to be larger than $B_{pop}$. If we calculate $S_{pop}$ for the two populations we have:

```{r}
Spec.barcelona <- mean(indvar.barcelona$B)/popvar.barcelona$B
Spec.barcelona
Spec.moia <- mean(indvar.moia$B)/popvar.moia$B
Spec.moia
```

Surprisingly, the degree of specialization in Moià seems slightly higher than that in Barcelona. To see whether this holds statistically, we can calculate the degree of specialization of each individual:
\[
S_i = \frac{B_i}{B_{pop}}
\]
which, in R, is:
```{r}
Spec.ind.barcelona <- indvar.barcelona$B/popvar.barcelona$B
Spec.ind.moia <- indvar.moia$B/popvar.moia$B
```

Finally, we compare this two vectors in a Wilcoxon test:
```{r}
wilcox.test(Spec.ind.barcelona, Spec.ind.moia)
```

Which tells us that those differences in individual specialization are not statistically significant.

### Measuring the degree of overlap between individuals
The idea of this section is to determine how much the niche of each individual overlaps with the niche of other individuals in the population. This can be done by calling function `nicheoverlap()` using `mode = "pairwise"`:
```{r}
O.barcelona <- nicheoverlap(diet.barcelona,D=dfood, mode="pairwise")
O.moia <- nicheoverlap(diet.moia,D=dfood, mode="pairwise")
```

These calls to `nicheoverlap()` return a symmetric square matrix with as many rows and columns as individuals in the resource use data frame. Each cell value in the symmetric matrix is the overlap between two individuals of the population. Using these matrices we can derive the average overlap in each population:
```{r}
mean(O.barcelona[lower.tri(O.barcelona)])
mean(O.moia[lower.tri(O.moia)])
```

We can also calculate the average overlap between each individual and the remaining individuals in its population:
```{r}
O.barcelona.ind <- (rowSums(O.barcelona)-1)/(nrow(O.barcelona)-1)
summary(O.barcelona.ind)
O.moia.ind <- (rowSums(O.moia)-1)/(nrow(O.moia)-1)
summary(O.moia.ind)
```

We substracted one in the numerator and denominator in order to exclude the target individual from the average (the overlap between a resource niche and itself is always one). Apparently, the individuals in Moià have a larger degree of overlap with individuals of their population than individuals in Barcelona. A non-parametric test seems to confirm this difference:
```{r}
wilcox.test(O.barcelona.ind, O.moia.ind)
```


## Bibliography