# 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 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")