Load packages and read in files

library(picante)
library(ape)
library(phylotools)
library(phytools)
library(terra)
library(geotax)
library(reshape2)
library(ggplot2)
library(gridExtra)

Look at tree and load files.

## Use vertlife tree for analysis ##
tree<- read.nexus("tree/output.nex")
consensus<-consensus.edges(tree, consensus.tree=consensus(tree,p=0.5))  #ape for majority rule consensus
plot(consensus)

dist<- read.csv("DistanceVertLife.csv", header=TRUE, row.names = 1) 
rpairs<- read.csv("RodPar_ChileFeb.csv", header=TRUE)
rhpairs<- read.csv("RodPar_ChileFeb_Hanta.csv", header=TRUE)
colnames(dist)<- row.names(dist)
incidence <- get_incidence_matrix(rpairs, returnDataFrame=FALSE) ## binary interaction matrix
incidence<- incidence[ ,colnames(dist)]
incidence<- subset(incidence, (rowSums(incidence)> 3)) # get rid of rodents with less than 3 ectoparasite associations
coef_all <- log_reg_boostrap(incidence, dist, 1000) # log reg coefs for all ectoparasites
coef_all
##     intercept    Std. Error       z value      Pr(>|z|)         2.5 % 
## -7.733556e-01  1.339613e-01 -5.772286e+00  4.835040e-06 -1.035915e+00 
##        97.5 %         slope    Std. Error       z value      Pr(>|z|) 
## -5.107962e-01 -2.120952e-02  2.463103e-03 -8.590565e+00  7.394558e-11 
##         2.5 %        97.5 % 
## -2.603712e-02 -1.638193e-02
coef <- sapply(1:nrow(incidence), function(x)
  log_reg_boostrap(incidence[x, ,drop = F], dist, 1000) ) #for each ectoparasite
colnames(coef)<- rownames(incidence) # change the v1,2,3 to the ectoparasite names
GD<- seq(0, 120,5)## scale for genetic distance
coefval <- t(coef)[ ,c(1,7)]
coef_allval<- c(coef_all[["intercept"]], coef_all[["slope"]] ) #subset the coefs to just slope/intercept
reg <- apply(coefval, 1, function(x) {
  prob_logit(x, GD)} ) # vector of probabilities for each ectoparasite
reg_all <- prob_logit(coef_allval, GD)

write.csv(coef_all, "coef_allpy-vert3.csv")
write.csv(coef, "coefpy-vert3.csv")
##reformat data to plot
plotreg<- cbind(reg, GD)
plotreg<-cbind(plotreg, reg_all)
plotreg<-as.data.frame(plotreg)
preg<-  melt(plotreg,  id.vars = 'GD', variable.name = 'series')

## do some weird stuff to get reg all line red
adj_names <- sort(setdiff(unique(preg$series), "reg_all" ))
gg_color_hue <- function(n) {
  hues = seq(0, 0, length = n)
  hcl(h = hues, l =80, c = 0)[1:n]
}
values <- gg_color_hue(length(adj_names))
names(values)<- adj_names
values<- c(values, c(reg_all="red"))

## plot 
ggphy<-ggplot(preg, aes(x=GD, y=value, colour=series))+
  geom_line(show.legend=F)+ theme_minimal()+
  scale_colour_manual(values=values)+
  xlab("Phylogenetic Distance between Rodents")+
  ylab("Probability of Sharing Ectoparasites")

ggphy