naturalis/trait-geo-diverse-ungulates

Comparative analysis of abiotic niches in Ungulates

The diversity of Ungulate abiotic niches

Project repository for the analysis of abiotic niches of Ungulates. Contains source code (R) and data files. The overall layout is as follows:

• script - reusable scripts in Perl and R
• data - data, including occurrences, GIS layers and shape files
• results - analysis output (e.g. maxent output)
• doc - (executable) documentation files, including RMarkdown

The source code that is being developed is intended to be portable such that analysis workflows can be prototyped on laptops, while longer running analyses can be executed in a cloud environment. In addition to the occurrences, other input data include GIS layers for bioclimatic variables, soil, topographic heterogeneity, and vegetation types. However, these are quite large files obtained from 3rd parties, and so we consider these immutable for this project. Hence, we may store some of these separately in a DropBox folder structure. From the combined input data types, we model niches (e.g. using maxent), the results of which we write to the results folder. These are then, subsequently, used for the following:

• to calculate niche overlap between pairs of species, such that we can then cluster species by ecological similarity. This is done in a variety of ways, including by computing projected niche overlap (Schoener's D) and by computing climatic 'trait' values (e.g. those from GIS pixels that coincide with raw occurrences, or GIS pixels in projected suitable habitats) and distances in trait space between species (Gower's D).
• to identify the niche dimensions that determine Ungulate distribution, i.e which GIS layers are the most important.

trait-geo-diverse-ungulates's Issues

Soil properties

In de literatuur kan ik nergens een species distribution model voor dieren vinden die soil properties meeneemt. De meeste modellen (met dieren) kijken naar klimaat en vegetatie. Vegetatie is natuurlijk een indirecte maat van de bodem kwaliteit. Ik denk dat voor het model wat we nu maken bodem niet belangrijk is maar als het model ook gebruikt wordt voor vegetatie dan zijn bodem kwaliteiten wel heel bepalend.

Data cleaning

De data zoals die uit de db kwamen moeten op de volgende wijze opgeschoond worden:

• alleen unieke records (dus niet dezelfde combinatie van lat/long meerdere keren)
• alleen records binnen de shape files van de species ranges (IUCN data)
• outliers moeten verwijderd worden, bijvoorbeeld door records die te veel afwijken van het "gemiddelde" eruit te gooien
• er moeten minimaal 10 records overblijven

Portable maken van deze repository

Om deze repository her-bruikbaar te maken voor verder onderzoek zou ik het volgende willen voorstellen:

• waarschijnlijk kunnen alle UNIEKE occurrences per soort in een CSV bestand, en dan al die bestanden in een map data in de root van deze repo
• vervolgens kunnen dan alle scripts daar met relatieve paden naar verwijzen (je hebt nu absolute paden in je code)

Niche modelleer vragen:

• wat is een goede benadering/package voor ensembles?
• wat is een redelijke manier om buffering te automatiseren (hoe veel % rondom occurrences)?
• wat is een redelijke manier om variable selection te automatiseren (#1)?
• hoe valideren we de modellen (AUC, congruentie van projectie en occurrences)? (#2)?
• hoe gaan we de data opschonen (#8)?
• hoe gaan we om met categorische variabelen (bodem, vegetatie)?

Variable selection

De set niet-gecorreleerde variabelen kan op verschillende manieren bepaald worden, bijvoorbeeld met de hand (te veel werk voor 220 soorten, lijkt me) of automatisch, door iteratief variabelen te verwijderen en AIC tests te doen. Het eindresultaat is dan dat er voor elke soort in principe een andere set variabelen geselecteerd zou kunnen worden. Is dat erg voor de vergelijkbaarheid? Bijvoorbeeld, als twee variabelen altijd met elkaar correleren (zeg, twee proxies die met temperatuur te maken hebben), dan zou het goed zijn voor de interpretatie dat we altijd dezelfde van de twee te selecteren. Maar misschien is het sowieso goed om ook modellen te bouwen waar alle soorten dezelfde set variabelen gebruiken? Punt van discussie...

5 minutes

hoe lang zou het duren om alles voor 5 arc minutes te draaien?

occurrence is met twee erren

Correlatie

Als we de correlatie tussen de abiotic lagen per species gaan berekenen weet ik niet of we de maxent resultaten wel met elkaar kunnen vergelijken (omdat je dan andere input in het model stopt).

downstream analysis SESHAT

http://seshatdatabank.info/data

README bestanden voor soorten

De meest elegante oplossing lijkt mij om een standaard template te maken, met daarin aangegeven wat de variabelen zijn. Bijvoorbeeld, waar de naam van de soort moet komen doen we: {{taxon_name}}. Dat is de standaard voor de template-taal mustache, waarvoor een R implementatie bestaat genaamd whisker.

nalopen modellen

De collectie van maxent modellen is wat onoverzichtelijk aan het worden. Kijkend in deze folder zijn er op dit moment de volgende modellen:

1. de modellen die @ElkeHendrix heeft gedraaid en die ik heb overgezet naar de nieuwe mappenstructuur. Deze modellen hebben allemaal de commit message 'adding raw occurrences and cached models'.
2. nieuwe modellen die ik heb gedraaid. Van deze modellen zijn er runs die zonder foutmeldingen zijn verlopen. Deze hebben het bericht 'recomputed with proper layer names'
3. een aantal runs hadden wel foutmeldingen, en daarvan is een veel voorkomende categorie de situatie waar een aantal input occurrences leiden tot NA predictor values, e.g. '1 (3.33%) of the presence points have NA predictor values'. Voor deze modellen moeten we de occurrences dus nalopen. Ik vermoed dat het komt doordat deze occurrences langs de kustlijn liggen en dus buiten het bereik van de layers komen.
4. voor een aantal kregen we andere foutmeldingen. Sus cebifrons en Sus philippensis zijn zwijnachtigen van de Filippijnen, en komen maar op bepaalde eilanden voor met een klein extent, dus de resampling werkt niet goed.
5. er zijn nog twee met onbegrijpelijke meldingen: 0c541db, 88d7151
6. een aantal modellen zijn 'invalid' volgens de validatiemethode

We moeten dus een aantal dingen nalopen:

• zijn alle modellen uit stap 1. inderdaad zonder foutmeldingen verlopen?
• handmatig de occurrences voor 3. en 4. checken
• uitvogelen wat er mis is voor 5.
• uitvogelen wat er mis is voor 6.

Checks:

• de volgorde van de namen van de layers (die in een vector stonden) klopte niet in relatie tot de volgorde waarin de lagen als TIF bestanden werden ingelezen. Vermoedelijk is het dus zo dat alle variable importance plots daardoor niet klopten.

  trait-geo-diverse-ungulates/doc/Niche_Overlap.Rmd

  Lines 103 to 131 in 2e8261f

  	files.names <- list.files(paste(REPO_HOME, "/data/GIS/5_deg", sep = ""))
  	# Turn the file names into layer names: strip the prefix, which might include
  	# the resolution, and strip the file extension
  	gis.layers.names <- files.names
  	gis.layers.names <- gsub('current_5arcmin_','',gis.layers.names)
  	gis.layers.names <- gsub('.tif','',gis.layers.names)
  	# Combine the layer names with those we've already read from BIOCLIM
  	gis.layers.names <- c(names(gis.layers),gis.layers.names)
  	# Iterate over files
  	for (i in 1:length(files.names)) {
  	# Stack with previously read layers
  	gis.layers <- stack(
  	gis.layers,
  	# Read as raster
  	raster(
  	# Construct file name
  	paste(REPO_HOME, "/data/GIS/5_deg/", files.names[i], sep = "")
  	)
  	)
  	}
  	# Apply all names
  	names(gis.layers) <- gis.layers.names

• klopt dit verhaal? Wat doen we hiermee?

  trait-geo-diverse-ungulates/doc/Niche_Overlap.Rmd

  Lines 257 to 263 in 2e8261f

  	Afterwards the occurrence dataset is split in k-fold partitions: a training
  	dataset containing 75% of the data and a test dataset containing 25% of the
  	data. The maxent model is constructed using the maxent function from the dismo
  	R package (Hijmans & Elith, 2013). The function extracts abiotic environmental
  	data for the training occurrence locations and 1000 random sampled background
  	locations, resulting in a model maxent object that can be used to predict which
  	other locations are suitable.

• dit moet anders. Als bepaalde soorten niet werken dan moeten ze structureel uit de lijst met taxa want op deze manier gaan we geheid ergens fouten maken met de volgorde

  trait-geo-diverse-ungulates/doc/Niche_Overlap.Rmd

  Line 295 in 2e8261f

  for (i in c(1:134, 136:150, 152:154)) {

• Vicagna vicugna klopt niet: de geslachtsnaam is incorrect en deze typefout heeft zich inmiddels door de rest van de bestanden en resultaten heen verspreid
• de syntax moet netter: we zouden middels styler (of anderszins) de code laten inspringen, dus niet zoals dit blok:

  trait-geo-diverse-ungulates/doc/Niche_Overlap.Rmd

  Lines 367 to 408 in 2e8261f

  	# combine two lists of valid and invalid models
  	output_AUC_valid <- data.frame(matrix(ncol = 3, nrow = length(list_species_model_high_accuracy)))
  	colnames(output_AUC_valid) <- c("taxon","trainingAUC","validation")
  	AUC.csv <- paste(REPO_HOME, "/results/maxent/AUCvalues.csv", sep="")
  	for (i in 1:length(list_species_model_high_accuracy)) {
  	open_species_model <- list_species_model_high_accuracy[[i]]
  	name <- names(list_species_model_high_accuracy[i])
  	name_underscore<- gsub( "_", " ", name)
  	trainingAUC<-open_species_model@results[[5,1]]
  	output_AUC_valid[i,] <- c(name_underscore, trainingAUC, "valid")
  	}
  	output_AUC_invalid <- data.frame(matrix(ncol = 3, nrow = length(list_species_model_low_accuracy)))
  	colnames(output_AUC_invalid) <- c("taxon","trainingAUC","validation")
  	for (i in 1:length(list_species_model_low_accuracy)) {
  	open_species_model <- list_species_model_low_accuracy[[i]]
  	name <- names(list_species_model_low_accuracy[i])
  	name_underscore<- gsub( "_", " ", name)
  	trainingAUC<-open_species_model@results[[5,1]]
  	output_AUC_invalid[i,] <- c(name_underscore, trainingAUC, "invalid")
  	}
  	combined_auc<- rbind(output_AUC_invalid, output_AUC_valid)
  	write.csv(combined_auc, file= AUC.csv)

In welke volgorde zijn de lagen van Totalstack.tif gestacked?

https://github.com/naturalis/trait-geo-diverse-ungulates/blob/master/script/Niche_Overlap.Rmd#L91-L101

model predictions werken nu niet meer?

In dit blok: https://github.com/naturalis/trait-geo-diverse-ungulates/blob/master/doc/1_maxent.rmd#L427-L492

Volgens mij kunnen we better dat hele clippen met die shape file achterwege laten want dan kunnen we tenminste nog Schoener's D berekenen. Nu crasht de berekening bij (nagenoeg) elke soort.

downstream analysis n-dimensional hypervolume

https://benjaminblonder.org/2014/02/26/the-ecological-niche-and-the-n-dimensional-hypervolume/

Bespreken donderdag

Ik heb een aantal punten die voor mij nog niet helemaal duidelijk zijn. Misschien kunnen we die donderdag bespreken:

• Hoe gaan we de OMI gebruiken --> nu gebruiken we alleen maar de tabel met de genormaliseerde waarden, maar mischien kunnen we de OMI ook nog plotten als contmap.

• Zullen we dit weekend meteen alles op 5 min runnen? dan zijn dat meteen de goede resultaten.

• In de results staat een clustering grafiek. Is die normale verdeling gemaakt over de volledige dataset of alleen over de domesticates? Of kan je in deze grafiek zien dat de APD voor domesticates kleiner (dus die rode streep) is dan die van de rest van de dataset.

• Gaan we ook een mantel test doen tussen onze boom en een echte phylogenetic trees?

• waar komen de phylogenetic trees vandaan?

• In het script niche traits bij het kopje r domestication haal je of de standardized.averages.csv op of je maakt er zelf een maar dat zijn de gemiddelden?

• In het script niche traits bij het kopje r phylo wordt een phylogenetic tree ingeladen en dan wordt er steeds een voorspellingsmodel gemaakt met verschillende abiotic variabelen? en het model met de laagste AIC is het beste in het voorspellen van domesticates gebaseerd op een x aantal variabelen? Ik snap hier eigenlijk ook niet waarom dit op de phylogenetic tree wordt gedaan en waarom niet gewoon op een dataframe met daarbij domesticate of wild?

• In het script niche traits bij het kopje r filter wordt gekeken of er nog andere niet gedomesticeerde soorten binnen de gedomesticeerde kenmerken vallen?

• afspraak maken met erica voor mid term gesprek en eindgesprek

falende gedomesticeerden

• Bos_taurus_primigenius 678e0d3
• Camelus_bactrianus 69f2336
• Equus_africanus

handwerk

Het valt me op dat de naam van de vicuna verkeerd is gespeld en dat er wat inconsistenties zijn in de kolommen van de gedomesticeerde data sets (quotes ja/nee, een eerste kolom met record numbers ja/nee). Ik vermoed dat dit het resultaat is van handwerk, wat dus niet goed reproduceerbaar is. Het lijkt me toch beter dat we dat niet doen: in principe moet elke stap opnieuw te doen zijn door scripts te runnen.

Scripts als rmarkdown

Bijvoorbeeld: https://github.com/naturalis/CourseComparativeMethods/blob/master/lecture1/Introduction_phylo.Rmd

coding style

Vanaf nu gaan we onze R code standaardizeren qua opmaak. Als volgt:

1. install.packages("styler")
2. RStudio > Tools > Addins > Browse addins...
3. selecteer 'styler ... pretty-print active file'

volledige namen van functies

De code in Niche_Overlap.Rmd laadt nu ongeveer twintig packages:

trait-geo-diverse-ungulates/script/Niche_Overlap.Rmd

Lines 16 to 37 in 2b3edbf

	library(raster, quietly = T)
	library(knitr, quietly = T)
	library(maxent, quietly = T)
	library(maps, quietly = T)
	library(rJava, quietly = T)
	library(maptools, quietly = T)
	library(jsonlite, quietly = T)
	library(caret, quietly = T)
	library(ENMeval, quietly = T)
	library(repmis, quietly = T)
	library(CoordinateCleaner, quietly = T)
	library(dismo, quietly = T)
	library(virtualspecies, quietly = T)
	library(sp, quietly = T)
	library(rgeos, quietly = T)
	library(ape, quietly = T)
	library(adehabitatMA, quietly = T)
	library(ade4, quietly = T)
	library(raster, quietly = T)
	library(SDMTools, quietly = T)
	library(factoextra, quietly = T)
	library(ecospat, quietly = T)

1. Zijn die allemaal nodig? Het is nu moeilijk te begrijpen wat we daadwerkelijk gebruiken.
2. Als ze inderdaad nodig zijn dan zou het handig zijn als we de volledige namen van de functies die we aanroepen in die packages gaan gebruiken.

Om een voorbeeld te geven: er zijn meerdere packages die een maxent functie hebben (dismo::maxent() en maxent::maxent()). Welke we nu aanroepen is afhankelijk van de volgorde waarin we packages laden, wat een nogal riskante vorm van action at a distance is.

Het is dus beter om expliciet aan te geven welke functie we aanroepen: dat maakt de code begrijpelijker en voorkomt bugs die puur komen door veranderingen in de volgorde waarin we packages laden.

Wilde voorouders

Veruit de meeste soorten zijn wild, dus de huidige waarnemingen komen in belangrijke mate overheen met hun niche gedurende het hele Holoceen. Voor gedomesticeerde soorten is dat niet zo: ze zijn gefokt voor verschillende landschapstypen en worden bijgevoerd en in stallen gehouden. Idealiter zouden we de modellen dus baseren op (subfossiele?) waarnemingen van hun wilde voorouders. Dat betekent misschien ook dat we dus paleoklimaat moeten meenemen?

Validatie

Hoe gaan we de kwaliteit van de modellen valideren? Er gaan vast situaties zijn waar de data gewoon niet goed genoeg zijn voor een bepaalde soort: te weinig waarnemingen, misidentificaties, "waarnemingen" in dierentuinen, verkeerde lat/lon coderingen, etc. Hoe gaan we objectief bepalen dat we een bepaalde soort echt niet mee kunnen nemen in de analyse?

manuscript

Link naar overleaf project: https://www.overleaf.com/project/5c7cfef8ac6a080f4fd4476a