<commit_before><img src="https://github.com/dsilvestro/PyRate/blob/master/pyrate_lib/PyRate_logo1024.png" align="left" width="80">  

# PyRate Tutorial \#1

#### Feb 2019

***
#### Contents
* [Generate PyRate input file (option 1)](https://github.com/dsilvestro/PyRate/blob/master/tutorials/pyrate_tutorial_1.md#generate-pyrate-input-file-option-1)  
* [Generate PyRate input file (option 2)](https://github.com/dsilvestro/PyRate/blob/master/tutorials/pyrate_tutorial_1.md#generate-pyrate-input-file-option-2)  
* [Estimation of speciation and extinction rates through time](https://github.com/dsilvestro/PyRate/blob/master/tutorials/pyrate_tutorial_1.md#estimation-of-speciation-and-extinction-rates-through-time)    

* [Return to Index](https://github.com/dsilvestro/PyRate/tree/master/tutorials#pyrate-tutorials---index) 
***

# Generate PyRate input file (option 1)
1. **Download fossil occurrences** for a clade from the [Paleobiology Database](https://paleobiodb.org/classic/displayDownloadGenerator). E.g. search for the genus *Canis* and save it as a cvs file, e.g. using the file name *Canis_pbdb_data.csv*. Before downloading the file, check the box "Show accepted names only" in the "Select by taxonomy" section, and uncheck the box "Include metadata at the beginning of the output" in the "Choose output options" section.

2. **Launch R** by opening a Terminal window and typing `R` then hit Enter. On Mac or Windows you can use the R GUI app or RStudio. 

3. **Load the *pyrate_utilities.r* file** e.g. by typing:  
  `source(".../PyRate-master/pyrate_utilities.r")`. Note that the full path of a file is here indicated as `.../` for simplicity. In most operating systems you can drag and drop the *pyrate_utilities.r* file onto the R console to paste the full path. 

4. **Check which species are extant today** (if any), as this information must be provided when running a diversification rate analysis. All species unless otherwise specified will be considered as extinct. Define a vector of extant species by typing in R:   
  `extant_dogs = c("Canis rufus","Canis lupus","Canis aureus","Canis latrans","Canis mesomelas","Canis anthus","Pseudalopex gymnocercus","Canis adustus","Canis familiaris")`

5. **Parse the raw data and generate PyRate input file**. Here we are going to use an automated function to extract fossil occurrence data from PBDB raw table and save it in a PyRate-compatible input file. Type in R:   
  `extract.ages.pbdb(file= ".../Canis_pbdb_data.csv",extant_species=extant_dogs)`  
  **IMPORTANT NOTE:** this function does not check for synonyms, typos, etc. You can specify two additional options.  The option `replicates` (by default set to 1) resamples the ages of each fossil occurrence from the respective temporal range. If you use `replicates=10`, ten replicated datasets will be saved in a single PyRate input file. These replicates can then be analyzed (individually) and the combined results will account for age uncertainties in the fossil record. The option `cutoff` can be used to remove occurrences with an age range greater than a given temporal range. For instance using `cutoff=10` will remove all occurrences in which the difference between max age and min age is larwger than 10 Myr. 

6. **Check PyRate input files.** The previous step generated two files: a text file (\*\_SpeciesList.txt) containing the list of all species included in the dataset and a python file (\*\_PyRate.py) with all occurrences formatted for a PyRate analysis. We can check the python file directly in PyRate to get a few summary statistics:  
  a. Open a Terminal window
  b. Browse to the Pyrate folder, e.g.:  
    `cd ".../PyRate-master"`
  c. Launch PyRate with the following arguments:  
  `python PyRate.py '.../Canis_pbdb_data_PyRate.py' -data_info`  
***

# Generate PyRate input file (option 2)
1. **Prepare fossil occurrence table.** You can prepare a table with fossil occurrence data in a text editor or a spreadsheet editor. The table must include 4 columns including 1) Taxon name, 2) Status specifying whether the taxon is "extinct" or "extant", 3) Minimum age, and 4) Maximum age. The table should have a header (first row) and **one row for each fossil occurrence** (see example below). Min and max ages represent the age ranges of each fossil occurrence, typically based on stratigraphic boundaries. One additional column can be included in the table indicating a trait value (e.g. estimated body mass) that can be used PyRate analyses.  


| Taxon\_name   | Status  | MinAge | MaxAge |
| ------------- |:-------------:| -----:| -----:|
Canis\_ferox | extinct | 4.9    | 10.3
Canis\_lepophagus | extinct | 4.9 | 10.3  
Canis\_ferox | extinct | 4.9 | 10.3  
Canis\_ferox | extinct | 4.9 | 10.3  
Canis\_ferox | extinct | 4.9 | 10.3  
Canis\_edwardii | extinct | 0.3 | 4.9  
Canis\_armbrusteri | extinct | 0.3 | 4.9  
Canis\_latrans | extant | 0.3 | 4.9  
Canis\_latrans | extant | 1.8 | 4.9  
Canis\_thooides | extinct | 1.8 | 4.9  

2. **Launch R** as explained above.  
3. **Load the *pyrate_utilities.r* file** as explained above.  
4. **Parse the raw data and generate PyRate input file**. Type in R:   
  `extract.ages(file="…/PyRate/example_files/Ursidae.txt", replicates=10)` This function includes `replicates` and `cutoff` options as the `extract.ages.pbdb()` function described above (option 1). 

### Accounting for age dependence in multiple occurrences from the same site
If a fossil assemblage or site contains several occurrences we should consider these occurrences as coeval even as we randomize their age within their temporal range. 
This can be done by specifying a "Site" column in the input data, with a number specifying the ID of the assemblage for each occurrence. For instance in the example below, the first three occurrences were found in the same site and their random age will be identical.


| Taxon\_name   | Status  | MinAge | MaxAge | Site
| ------------- |:-------------:| -----:| -----:| -----:|
Canis\_ferox | extinct | 4.9    | 10.3 | 1
Canis\_lepophagus | extinct | 4.9 | 10.3  | 1
Canis\_ferox | extinct | 4.9 | 10.3  | 1
Canis\_ferox | extinct | 4.9 | 10.3  | 2
Canis\_ferox | extinct | 4.9 | 5.3  | 3
Canis\_edwardii | extinct | 0.3 | 4.9  | 4
Canis\_armbrusteri | extinct | 0.3 | 4.9 | 4 
Canis\_latrans | extant | 0.3 | 4.9  | 4 
Canis\_latrans | extant | 1.8 | 4.9  | 5
Canis\_thooides | extinct | 1.8 | 4.9  | 5

*pyrate_utilities.r* (`extract.ages()`) will automatically resample fossil ages by site when a column named "Site" is included in the input table. **We thank <NAME> and <NAME> for input on this!**

***
# Estimation of speciation and extinction rates through time

## Defining the preservation model
**Non-homogeneous Poisson process of preservation (NHPP)**. By default, PyRate assumes an NHPP model in which preservation rates change during the lifespan of each lineage following a bell-shaped distribution (Silvestro et al. 2014 Syst Biol). 

**Homogeneous Poisson process (HPP)**. Alternatively, an HPP model, in which the preservation rate is constant through time, is also available, using the command `-mHPP`:

`python PyRate.py .../Canis_pbdb_data_PyRate.py -mHPP`

**Time-variable Poisson process (TPP)**. An alternative model of preservation assumes that preservation rates are constant within a predefined time frame, but can vary across time frames (e.g. geological epochs). This model is particularly useful if we expect rate heterogeneity to occur mostly through time, rather than among lineages.

We can set up a model in which preservation rates are estimated independently within each geological epoch by using the command `-qShift` and providing a file containing the times that delimit the different epochs (an example is provided in the example file _PyRate-master/example\_files/epochs\_q.txt_):  

`python PyRate.py .../Canis_pbdb_data_PyRate.py -qShift .../epochs_q.txt`  

The default prior on the vector of preservation rates is a single gamma distribution with shape = 1.5 and rate = 1.5. These parameter values can be changed using the command `-pQ`, for example `-pQ 2 0.1` changes the shape to 2 and the rate to 0.1.
Alternatively, the rate parameter of the prior can be estimated from the data to reduce the subjectivity of the choice and allow the prior to better adapt to the data. This option is available by setting the rate parameter to 0, e.g.: 

`python PyRate.py .../Canis_pbdb_data_PyRate.py -qShift .../epochs_q.txt -pQ 1.5 0`

When the rate is set to zero, PyRate assigns a vague exponential hyper-prior to the rate and samples the rate along with all other model parameters. This approach is recommended whenever multiple preservation rates are estimated


**Gamma model of rate heterogeneity**. NHPP, HPP, and TPP models can all be coupled with a Gamma model of rate heterogeneity, which enables us to account for heterogeneity in the preservation rate across lineages. This option only adds a single parameter to the model and should be used for all empirical data sets. To set the Gamma model we add the flag `-mG`:

`python PyRate.py .../Canis_pbdb_data_PyRate.py -mG` [NHPP model]  
`python PyRate.py .../Canis_pbdb_data_PyRate.py -mHPP -mG` [HPP model]
`python PyRate.py .../Canis_pbdb_data_PyRate.py -qShift .../epochs_q.txt -mG`  


### Model testing across preservation models
A maximum likelihood test across (described [here](https://www.biorxiv.org/content/10.1101/316992v1)) is available to assess which of NHPP, HPP, or TPP is best supported by the data. To run the test you need to provide the input file and the file providing the times of rate shift for the TPP model and add the flag `-PPmodeltest`

``python PyRate.py .../Canis_pbdb_data_PyRate.py -qShift .../epochs_q.txt -PPmodeltest``

Note that the Gamma model is not tested here, but we recommend to add it to the best model among NHPP, HPP, or TPP as selected by the model testing. 

## Analysis setup
Here we describe the main settings for a standard analysis of fossil occurrence data using PyRate. The analysis will estimate:

1. **origination and extinction times** of each lineage  <commit_msg>fixed documentation hyper-prior on preservation rates<commit_after><img src="https://github.com/dsilvestro/PyRate/blob/master/pyrate_lib/PyRate_logo1024.png" align="left" width="80">  

# PyRate Tutorial \#1

#### Feb 2019

***
#### Contents
* [Generate PyRate input file (option 1)](https://github.com/dsilvestro/PyRate/blob/master/tutorials/pyrate_tutorial_1.md#generate-pyrate-input-file-option-1)  
* [Generate PyRate input file (option 2)](https://github.com/dsilvestro/PyRate/blob/master/tutorials/pyrate_tutorial_1.md#generate-pyrate-input-file-option-2)  
* [Estimation of speciation and extinction rates through time](https://github.com/dsilvestro/PyRate/blob/master/tutorials/pyrate_tutorial_1.md#estimation-of-speciation-and-extinction-rates-through-time)    

* [Return to Index](https://github.com/dsilvestro/PyRate/tree/master/tutorials#pyrate-tutorials---index) 
***

# Generate PyRate input file (option 1)
1. **Download fossil occurrences** for a clade from the [Paleobiology Database](https://paleobiodb.org/classic/displayDownloadGenerator). E.g. search for the genus *Canis* and save it as a cvs file, e.g. using the file name *Canis_pbdb_data.csv*. Before downloading the file, check the box "Show accepted names only" in the "Select by taxonomy" section, and uncheck the box "Include metadata at the beginning of the output" in the "Choose output options" section.

2. **Launch R** by opening a Terminal window and typing `R` then hit Enter. On Mac or Windows you can use the R GUI app or RStudio. 

3. **Load the *pyrate_utilities.r* file** e.g. by typing:  
  `source(".../PyRate-master/pyrate_utilities.r")`. Note that the full path of a file is here indicated as `.../` for simplicity. In most operating systems you can drag and drop the *pyrate_utilities.r* file onto the R console to paste the full path. 

4. **Check which species are extant today** (if any), as this information must be provided when running a diversification rate analysis. All species unless otherwise specified will be considered as extinct. Define a vector of extant species by typing in R:   
  `extant_dogs = c("Canis rufus","Canis lupus","Canis aureus","Canis latrans","Canis mesomelas","Canis anthus","Pseudalopex gymnocercus","Canis adustus","Canis familiaris")`

5. **Parse the raw data and generate PyRate input file**. Here we are going to use an automated function to extract fossil occurrence data from PBDB raw table and save it in a PyRate-compatible input file. Type in R:   
  `extract.ages.pbdb(file= ".../Canis_pbdb_data.csv",extant_species=extant_dogs)`  
  **IMPORTANT NOTE:** this function does not check for synonyms, typos, etc. You can specify two additional options.  The option `replicates` (by default set to 1) resamples the ages of each fossil occurrence from the respective temporal range. If you use `replicates=10`, ten replicated datasets will be saved in a single PyRate input file. These replicates can then be analyzed (individually) and the combined results will account for age uncertainties in the fossil record. The option `cutoff` can be used to remove occurrences with an age range greater than a given temporal range. For instance using `cutoff=10` will remove all occurrences in which the difference between max age and min age is larwger than 10 Myr. 

6. **Check PyRate input files.** The previous step generated two files: a text file (\*\_SpeciesList.txt) containing the list of all species included in the dataset and a python file (\*\_PyRate.py) with all occurrences formatted for a PyRate analysis. We can check the python file directly in PyRate to get a few summary statistics:  
  a. Open a Terminal window
  b. Browse to the Pyrate folder, e.g.:  
    `cd ".../PyRate-master"`
  c. Launch PyRate with the following arguments:  
  `python PyRate.py '.../Canis_pbdb_data_PyRate.py' -data_info`  
***

# Generate PyRate input file (option 2)
1. **Prepare fossil occurrence table.** You can prepare a table with fossil occurrence data in a text editor or a spreadsheet editor. The table must include 4 columns including 1) Taxon name, 2) Status specifying whether the taxon is "extinct" or "extant", 3) Minimum age, and 4) Maximum age. The table should have a header (first row) and **one row for each fossil occurrence** (see example below). Min and max ages represent the age ranges of each fossil occurrence, typically based on stratigraphic boundaries. One additional column can be included in the table indicating a trait value (e.g. estimated body mass) that can be used PyRate analyses.  


| Taxon\_name   | Status  | MinAge | MaxAge |
| ------------- |:-------------:| -----:| -----:|
Canis\_ferox | extinct | 4.9    | 10.3
Canis\_lepophagus | extinct | 4.9 | 10.3  
Canis\_ferox | extinct | 4.9 | 10.3  
Canis\_ferox | extinct | 4.9 | 10.3  
Canis\_ferox | extinct | 4.9 | 10.3  
Canis\_edwardii | extinct | 0.3 | 4.9  
Canis\_armbrusteri | extinct | 0.3 | 4.9  
Canis\_latrans | extant | 0.3 | 4.9  
Canis\_latrans | extant | 1.8 | 4.9  
Canis\_thooides | extinct | 1.8 | 4.9  

2. **Launch R** as explained above.  
3. **Load the *pyrate_utilities.r* file** as explained above.  
4. **Parse the raw data and generate PyRate input file**. Type in R:   
  `extract.ages(file="…/PyRate/example_files/Ursidae.txt", replicates=10)` This function includes `replicates` and `cutoff` options as the `extract.ages.pbdb()` function described above (option 1). 

### Accounting for age dependence in multiple occurrences from the same site
If a fossil assemblage or site contains several occurrences we should consider these occurrences as coeval even as we randomize their age within their temporal range. 
This can be done by specifying a "Site" column in the input data, with a number specifying the ID of the assemblage for each occurrence. For instance in the example below, the first three occurrences were found in the same site and their random age will be identical.


| Taxon\_name   | Status  | MinAge | MaxAge | Site
| ------------- |:-------------:| -----:| -----:| -----:|
Canis\_ferox | extinct | 4.9    | 10.3 | 1
Canis\_lepophagus | extinct | 4.9 | 10.3  | 1
Canis\_ferox | extinct | 4.9 | 10.3  | 1
Canis\_ferox | extinct | 4.9 | 10.3  | 2
Canis\_ferox | extinct | 4.9 | 5.3  | 3
Canis\_edwardii | extinct | 0.3 | 4.9  | 4
Canis\_armbrusteri | extinct | 0.3 | 4.9 | 4 
Canis\_latrans | extant | 0.3 | 4.9  | 4 
Canis\_latrans | extant | 1.8 | 4.9  | 5
Canis\_thooides | extinct | 1.8 | 4.9  | 5

*pyrate_utilities.r* (`extract.ages()`) will automatically resample fossil ages by site when a column named "Site" is included in the input table. **We thank <NAME> and <NAME> for input on this!**

***
# Estimation of speciation and extinction rates through time

## Defining the preservation model
**Non-homogeneous Poisson process of preservation (NHPP)**. By default, PyRate assumes an NHPP model in which preservation rates change during the lifespan of each lineage following a bell-shaped distribution (Silvestro et al. 2014 Syst Biol). 

**Homogeneous Poisson process (HPP)**. Alternatively, an HPP model, in which the preservation rate is constant through time, is also available, using the command `-mHPP`:

`python PyRate.py .../Canis_pbdb_data_PyRate.py -mHPP`

**Time-variable Poisson process (TPP)**. An alternative model of preservation assumes that preservation rates are constant within a predefined time frame, but can vary across time frames (e.g. geological epochs). This model is particularly useful if we expect rate heterogeneity to occur mostly through time, rather than among lineages.

We can set up a model in which preservation rates are estimated independently within each geological epoch by using the command `-qShift` and providing a file containing the times that delimit the different epochs (an example is provided in the example file _PyRate-master/example\_files/epochs\_q.txt_):  

`python PyRate.py .../Canis_pbdb_data_PyRate.py -qShift .../epochs_q.txt`  

The default prior on the vector of preservation rates is a single gamma distribution with shape = 1.5 and rate = 1.5. These parameter values can be changed using the command `-pP`, for example `-pP 2 0.1` changes the shape to 2 and the rate to 0.1.
Alternatively, the rate parameter of the prior can be estimated from the data to reduce the subjectivity of the choice and allow the prior to better adapt to the data. This option is available by setting the rate parameter to 0, e.g.: 

`python PyRate.py .../Canis_pbdb_data_PyRate.py -qShift .../epochs_q.txt -pP 1.5 0`

When the rate is set to zero, PyRate assigns a vague exponential hyper-prior to the rate and samples the rate along with all other model parameters. This approach is recommended whenever multiple preservation rates are estimated


**Gamma model of rate heterogeneity**. NHPP, HPP, and TPP models can all be coupled with a Gamma model of rate heterogeneity, which enables us to account for heterogeneity in the preservation rate across lineages. This option only adds a single parameter to the model and should be used for all empirical data sets. To set the Gamma model we add the flag `-mG`:

`python PyRate.py .../Canis_pbdb_data_PyRate.py -mG` [NHPP model]  
`python PyRate.py .../Canis_pbdb_data_PyRate.py -mHPP -mG` [HPP model]
`python PyRate.py .../Canis_pbdb_data_PyRate.py -qShift .../epochs_q.txt -mG`  


### Model testing across preservation models
A maximum likelihood test across (described [here](https://www.biorxiv.org/content/10.1101/316992v1)) is available to assess which of NHPP, HPP, or TPP is best supported by the data. To run the test you need to provide the input file and the file providing the times of rate shift for the TPP model and add the flag `-PPmodeltest`

``python PyRate.py .../Canis_pbdb_data_PyRate.py -qShift .../epochs_q.txt -PPmodeltest``

Note that the Gamma model is not tested here, but we recommend to add it to the best model among NHPP, HPP, or TPP as selected by the model testing. 

## Analysis setup
Here we describe the main settings for a standard analysis of fossil occurrence data using PyRate. The analysis will estimate:

1. **origination and extinction times** of each lineage  
{'fraction_non_alphanumeric': 0.09444280224652675, 'fraction_numerical': 0.017341610010838508, 'mean_word_length': 3.9401314188367, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 14, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1054471', 'n_tokens_mistral': 6362, 'n_tokens_neox': 5769, 'n_words': 2679}