RPhylopars

# Load libraries
library(phytools)

## Loading required package: ape

## Loading required package: maps

library(Rphylopars)

## You may need to install phylolm before Rphylopars will work. If so, uncomment and run this.
# devtools::install_github("lamho86/phylolm")

Run Phylopars

Run for each month + tree file combination across all sites.

# set working directory. This will need to be updated
setwd('C:/Users/april/OneDrive/Desktop/Covid/phylopars_updated_7-22-21/netherlands')

# read tree file
tree <- read.tree(file = "./RAxML_TreeTime_non_recombinant.newick")

# read data files and run phylopars
trait_data <- read.csv("jan.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
jan <- PPE$anc_recon 

trait_data <- read.csv("feb.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
feb <- PPE$anc_recon 

trait_data <- read.csv("mar.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
mar <- PPE$anc_recon 

trait_data <- read.csv("apr.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
apr <- PPE$anc_recon 

trait_data <- read.csv("may.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
may <- PPE$anc_recon 

trait_data <- read.csv("jun.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
jun <- PPE$anc_recon 

trait_data <- read.csv("jul.csv",h=T) 
PPE_jul <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
jul <- PPE$anc_recon 

trait_data <- read.csv("aug.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
aug <- PPE$anc_recon 

trait_data <- read.csv("sep.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
sep <- PPE$anc_recon 

trait_data <- read.csv("oct.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
oct <- PPE$anc_recon 

trait_data <- read.csv("nov.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
nov <- PPE$anc_recon 

trait_data <- read.csv("dec.csv",h=T) 
PPE <- phylopars(trait_data,tree,model="BM",phylo_correlated = TRUE, pheno_correlated = TRUE) 
dec <- PPE$anc_recon 

# combine and write to csv
comb <- rbind(jan, feb, mar, apr, may, jun, jul, aug, sep, oct, nov, dec)
write.csv(comb, file = "RAxML_TreeTime_non_recombinant.csv")