install.packages('ape', repos='http://cran.wustl.edu/')
install.packages('seqinr', repos='http://cran.wustl.edu/')
library(ape)
library(seqinr)

data <- read.table('Salmon_data.txt', header=TRUE)
newick <- read.tree('Salmon.nwk')
x<- data$body_length
names(x)<- data$species
y<-data$Egg_weight
names(y)<- data$species
original_correlation = cor(pic(x, newick), pic(y, newick))
original_correlation

## Code for including phylogenetic uncertainty when the uncertainty
## is huge - as in the phylogeny is completely unknown
B=1000
n = nrow(data)
rcorrelation = array(0, dim=c(B))
z = array(0, dim=c(B))
for(i in 1:B){
rnewick <-  rtree(n, rooted = TRUE, tip.label = data$species, br = runif)  
rcorrelation[i] = cor(pic(x, rnewick), pic(y, rnewick))
z[i] = .5*log((1 + rcorrelation[i])/(1 - rcorrelation[i]))
}
average_z = mean(z)
average_correlation = tanh(average_z)
var_z = (1/(n-4)) + var(z)
z_l = average_z - 1.96*sqrt(var_z)
z_u = average_z + 1.96*sqrt(var_z)
r_l = tanh(z_l)
r_u = tanh(z_u)
c(r_l, average_correlation, r_u)