---
title: "Case studies (GEB_739_sm_AppendixS2-S5)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Case studies (GEB_739_sm_AppendixS2-S5)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, echo=FALSE}
library(knitr)
opts_chunk$set(
  warning = FALSE, 
  message = FALSE, 
  fig.width = 4,   
  fig.height = 5,   
  collapse = TRUE,
  comment = "#>"
  )
```

Supporting information in
<a href="https://onlinelibrary.wiley.com/doi/10.1111/j.1466-8238.2011.00739.x" target="_blank">
Valcu, M., Dale, J., and Kempenaers, B. (2012). 
rangeMapper: a platform for the study of macroecology of life-history traits.
Global Ecology and Biogeography 21, 945-951.
</a>


###  General project setup

We setup a project using a hexagonal canvas with a cell size of 500 km.
The project is set `in-memory` but for a real case study you would like to set `path` to a persistent location on disk.  
We'll use the `wrens` dataset which is part of the package. 


```{r}
require(rangeMapper)
require(sf)
require(data.table)
require(glue)
require(ggplot2)
require(viridis)

wrens = read_wrens()
wrens$breeding_range_area = st_area(wrens)

con = rmap_connect()

rmap_add_ranges(con, x = wrens, ID = 'sci_name')
rmap_prepare(con, 'hex', cellsize = 500)
rmap_add_bio(con, wrens, 'sci_name')

```

### Raw data: wrens breeding range distribution and life history


```{r}

ggplot() + 
  geom_sf(data = rmap_to_sf(con, 'wkt_canvas') , color = 'grey80', fill = NA) +
  geom_sf(data = wrens, fill = NA) 

head(wrens,3)

```


## Case study 1: Different biodiversity hotspots and their congruence.

We describe biodiversity hotspots based on of three  avian diversity parameters: **total species richness**, **endemic species richness** and **relative body mass diversity**. 


### 1. Set parameters

```{r}
P_richness  = 0.75 # species richness quantile
P_bodymass  = 0.50 # CV body mass quantile
P_endemics  = 0.35 # endemic species richness quantile

```

### 2. Primary maps and thresholds


```{r}
rmap_save_map(con)  
# rmap_save_map with no arguments other than `con` saves a species_richness map.

CV_Mass <- function(x) (sd(log(x),na.rm = TRUE)/mean(log(x),na.rm = TRUE))
rmap_save_map(con, fun = CV_Mass, src='wrens',v = 'body_mass', dst='CV_Mass')

```

Thresholds are computed using the parameters defined in **1** and the maps saved at **2**. 

```{r}
sr = rmap_to_sf(con, "species_richness")
sr_threshold = quantile(sr$species_richness, probs = P_richness, na.rm = TRUE)

es_threshold = quantile(wrens$breeding_range_area, probs = P_endemics, na.rm = TRUE)

bmr = rmap_to_sf(con, "CV_Mass")
bmr_threshold = quantile(bmr$V1_body_mass, probs = P_bodymass, na.rm = TRUE)

```

### 3. Congruence subsets and congruence maps

#### 3.1 Subsets

```{r}
rmap_save_subset(con,'sr_threshold', species_richness = paste('species_richness     >', sr_threshold) )
rmap_save_subset(con,'es_threshold', wrens            = paste('breeding_range_area <=', es_threshold) )
rmap_save_subset(con,'bmr_threshold', CV_Mass         = paste('V1_body_mass        >=', bmr_threshold))

rmap_save_subset(con, "cumul_congruence_threshold",
    species_richness = paste('species_richness >', sr_threshold),
    wrens            = paste('breeding_range_area <=', es_threshold), 
    CV_Mass          = paste('body_mass >=', bmr_threshold) 
    )
```

#### 3.2 Threshold Maps

```{r}
rmap_save_map(con, subset = 'sr_threshold', dst = 'Species_richness_hotspots')
rmap_save_map(con, subset = 'es_threshold', dst = 'Endemics_hotspots')
rmap_save_map(con, subset = 'bmr_threshold', dst = 'Body_mass_diversity_hotspots')
rmap_save_map(con, subset = 'cumul_congruence_threshold', dst = 'Cumul_congruence_hotspots')


```

### 4. Maps: load and display


```{r, fig.height = 10, fig.width = 8}
study_area = rmap_to_sf(con, 'species_richness')  %>% st_union
bmr = rmap_to_sf(con, pattern = 'hotspots')  %>% 
      melt(id.vars = c('geometry', 'cell_id') )  %>% 
      st_as_sf
bmr$variable = bmr$variable %>% gsub('species_richness_|_hotspots', '', .)      

ggplot() + 
  facet_wrap(~variable) + 
  geom_sf(data = study_area ) + 
  geom_sf(data = bmr, aes(fill = value),  size= 0.05) + 
  scale_fill_gradientn(colours = viridis(10, option = 'E'), na.value= 'grey80') + 
  guides(fill=guide_legend(title='Wren\nspecies')) +
  ggtitle("Hotspots") +
  theme_bw()


```


## Case study 2: Geographical variation of the range size`~` body size slope

### 1. The range size`~` body size slope map


```{r}

lm_slope = function (x) {
  lm(scale(log(breeding_range_area)) ~ scale(male_tarsus), x)  %>% 
  summary %>% coefficients %>% data.frame %>% .[-1, ] 
  }


rmap_save_map(con, fun = lm_slope, src='wrens', dst='slope_area_body_mass')

```

### 2. Maps: load and display
```{r}
m = rmap_to_sf(con, 'slope_area_body_mass')

ggplot(m) + 
  geom_sf(aes(fill = Estimate),  size= 0.05, show.legend = TRUE) + 
  scale_fill_gradientn(colours = viridis(10, option = 'E') ) + 
  ggtitle("Range size ~ Body size slope")

```


## Case study 3: The influence of cell size on body size `~` species richness slope

### 1. assemblage level median body size `~` species richness slope for varying cell sizes. 


```{r}

cellSizes = seq(from = 700, to = 1500, length.out = 5)

FUN = function(g) {
  options(rmap.verbose = FALSE)
  
  con = rmap_connect()
  rmap_add_ranges(con, x = wrens, ID = 'sci_name')
  rmap_prepare(con, 'hex', cellsize=g)
  rmap_add_bio(con, wrens, 'sci_name')
  rmap_save_map(con)
  rmap_save_map(con, fun = 'median', src='wrens', v = 'male_tarsus', dst='median_male_tarsus')
  m = rmap_to_sf(con)

  # lm at assemblage level
  o = lm(scale(log(median_male_tarsus)) ~ sqrt(species_richness), m)  %>% 
        summary %>% coefficients %>% data.frame %>% .[-1, ]

  o$cell_size = g

  options(rmap.verbose = TRUE)

  o

  }


o = lapply(cellSizes, FUN)   %>% rbindlist

```

### 2. Plot regression parameters for different cell sizes

Most of the variation here is due to spatial autocorrelation, a proper analysis requires a spatial model. 


```{r}

ggplot(o, aes(x = cell_size, y = Estimate)) +
    geom_point() +
    theme_bw()

```


## Case study 4: The influence of range size on the relationship between species richness and body size 

### 1. Set parameters

```{r}
quant = seq(0.05, 1, 0.1) 
Q = quantile(wrens$breeding_range_area,  probs = quant)
range_classes = data.frame(area = Q, quant =  quant )
W = 4   # size of the moving window
maxn = nrow(range_classes) - W

range_classes

```
### 2. Make `subset` tables for multiple range size intervals. 

```{r}

subsets = paste0('area_subset_',1:maxn )

for(i in 1:maxn ) {
  rmap_save_subset(con, subsets[i], 
    wrens = glue("breeding_range_area BETWEEN 
            {range_classes[i,     'area']} AND 
            {range_classes[i+W, 'area']  }") )
  }

```

### 2. Make median body mass `maps` for all `subset` tables. 


```{r}

maps = paste0('body_size_', subsets)

for(i in 1:maxn ) {
  rmap_save_map(con, subset = subsets[i], dst = maps[i],  
                fun = 'median', src='wrens', v = 'male_tarsus')
    }

```

### 3. Get `maps` and run Run `assemblage body size ~ species_richness` regression


```{r}

m = rmap_to_sf(con, pattern = 'richness|body')  %>% setDT
m = melt(m, measure.vars = patterns("median") )

x = m[, { 
      
      fm = lm( log10(value) ~  sqrt(species_richness) )
      data.table( Estimate = coefficients(fm)[2], t( confint(fm)[2, ]) )
      
      } , by = variable]

x[, Quantiles :=  quant[1:maxn]  ]

```

### 4. Plot regression slope for different quantile-based range size intervals


```{r}
ggplot(x, aes(x = Quantiles, y = Estimate) ) +
    geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), width= 0) +
    geom_line() +
    geom_point() +
    theme_bw()

```