# County PD, MPD, MNND, ses.PD, ses.MPD, ses.MNND #
library(PhyloMeasures)
library(ape)
wd <- "/species_number/"
setwd(wd)

co <- read.csv("american_matrix_MASTER.csv") #everything in 48 states.
row.names(co) <- co$X
co <- co[, -1]

setwd("trees")

# multiple trees:
trees <- read.tree("Tallen2015_ALL.tre") # Best ML tree + 100 bootstrap replicates, all dated.
tree <- trees[[1]]
trees <- lapply(trees, drop.tip, tip = tree$tip.label[!tree$tip.label %in% colnames(co)]) # drop taxa not present in counties

co2 <- co[, colnames(co)[colnames(co) %in% tree$tip.label]] #subset data for faster processing: remove taxa not on tree
co2 <- co2[, order(names(co2))] #just in case

ontree <- rowSums(co2) 
total <- rowSums(co)

# get native status
setwd(wd)
status <- read.csv("sp2.csv")
ss.n <- as.character(status$sp[status$old_status == "N" | status$old_status == "N?"]) # all native
ss.i <- as.character(status$sp[status$old_status == "I" | status$old_status == "I?"]) # all introduced
ss.ni <- ss.i[!ss.i %in% status$sp[status$old_invasive == "Y"]] # all non-invasive introduced
ss.ii <- ss.i[ss.i %in% status$sp[status$old_invasive == "Y"]] # all invasive introduced

# get total species per category
total.n <- rowSums(co[, colnames(co) %in% ss.n]) # all native
total.i <- rowSums(co[, colnames(co) %in% ss.i]) # all introduced
total.ni <- rowSums(co[, colnames(co) %in% ss.ni]) # all non-invasive introduced
total.ii <- rowSums(co[, colnames(co) %in% ss.ii]) # all invasive introduced

co2.n <- co2[, colnames(co2) %in% ss.n]
co2.i <- co2[, colnames(co2) %in% ss.i]
co2.ni <- co2[, colnames(co2) %in% ss.ni]
co2.ii <- co2[, colnames(co2) %in% ss.ii]

ontree.n <- rowSums(co2.n)
ontree.i <- rowSums(co2.i)
ontree.ni <- rowSums(co2.ni)
ontree.ii <- rowSums(co2.ii)

# calculations using PhyloMeasures:
require(doSNOW)
require(doParallel)
cl <- makeCluster(8) #careful with memory requirements
registerDoParallel(cl)

setwd(wd)
setwd("output_alpha")

foreach(i = 1:length(trees),.packages=c('PhyloMeasures', 'ape')) %dopar% {
  setwd(wd)
  setwd("output_alpha")
  tree <- trees[[i]]
  # PD
  PD <- pd.query(tree, co2, standardize = FALSE)
  PD.s <- pd.query(tree, co2, standardize = TRUE)
  PD.p <- pd.pvalues(tree, co2)
  # MPD
  MPD <- mpd.query(tree, co2, standardize = FALSE)
  MPD.s <- mpd.query(tree, co2, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree, co2)
  # MNTD
  MNTD <- mntd.query(tree, co2, standardize = FALSE)
  MNTD.s <- mntd.query(tree, co2, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree, co2)

  out <- cbind(fips = row.names(co2), total, ontree, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("county_diversity.", i-1, ".csv"), row.names = FALSE)
  
  # by status:
  tree.n <- drop.tip(tree, colnames(co2)[!colnames(co2) %in% ss.n])
  tree.i <- drop.tip(tree, colnames(co2)[!colnames(co2) %in% ss.i])
  tree.ni <- drop.tip(tree, colnames(co2)[!colnames(co2) %in% ss.ni])
  tree.ii <- drop.tip(tree, colnames(co2)[!colnames(co2) %in% ss.ii])

  # natives:
  # 1. null = all taxa in US regardless of status
  co2.n <- co2
  co2.n[, colnames(co2.n)[!colnames(co2.n) %in% ss.n]] <- 0
  
  # PD
  PD <- pd.query(tree, co2.n, standardize = FALSE)
  PD.s <- pd.query(tree, co2.n, standardize = TRUE)
  PD.p <- pd.pvalues(tree, co2.n)
  # MPD
  MPD <- mpd.query(tree, co2.n, standardize = FALSE)
  MPD.s <- mpd.query(tree, co2.n, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree, co2.n)
  # MNTD
  MNTD <- mntd.query(tree, co2.n, standardize = FALSE)
  MNTD.s <- mntd.query(tree, co2.n, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree, co2.n)

  out <- cbind(fips = row.names(co2.n), total.n, ontree.n, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("NATIVE.county_diversity_ALLnull.", i-1, ".csv"), row.names = FALSE)
  
  # 2. null = all NATIVE taxa in US
  co2.n <- co2[, colnames(co2) %in% ss.n]

  # PD
  PD <- pd.query(tree.n, co2.n, standardize = FALSE)
  PD.s <- pd.query(tree.n, co2.n, standardize = TRUE)
  PD.p <- pd.pvalues(tree.n, co2.n)
  # MPD
  MPD <- mpd.query(tree.n, co2.n, standardize = FALSE)
  MPD.s <- mpd.query(tree.n, co2.n, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree.n, co2.n)
  # MNTD
  MNTD <- mntd.query(tree.n, co2.n, standardize = FALSE)
  MNTD.s <- mntd.query(tree.n, co2.n, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree.n, co2.n)

  out <- cbind(fips = row.names(co2.n), total.n, ontree.n, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("NATIVE.county_diversity_NATnull.", i-1, ".csv"), row.names = FALSE)
 
  # introduced:
  # 1. null = all taxa in US regardless of status
  co2.i <- co2
  co2.i[, colnames(co2.i)[!colnames(co2.i) %in% ss.i]] <- 0

  # PD
  PD <- pd.query(tree, co2.i, standardize = FALSE)
  PD.s <- pd.query(tree, co2.i, standardize = TRUE)
  PD.p <- pd.pvalues(tree, co2.i)
  # MPD
  MPD <- mpd.query(tree, co2.i, standardize = FALSE)
  MPD.s <- mpd.query(tree, co2.i, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree, co2.i)
  # MNTD
  MNTD <- mntd.query(tree, co2.i, standardize = FALSE)
  MNTD.s <- mntd.query(tree, co2.i, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree, co2.i)

  out <- cbind(fips = row.names(co2.i), total.i, ontree.i, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("INTRO.county_diversity_ALLnull.", i-1, ".csv"), row.names = FALSE)
  
  # 2. null = all INTRODUCED taxa in US
  co2.i <- co2[, colnames(co2) %in% ss.i]  
  
  # PD
  PD <- pd.query(tree.i, co2.i, standardize = FALSE)
  PD.s <- pd.query(tree.i, co2.i, standardize = TRUE)
  PD.p <- pd.pvalues(tree.i, co2.i)
  # MPD
  MPD <- mpd.query(tree.i, co2.i, standardize = FALSE)
  MPD.s <- mpd.query(tree.i, co2.i, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree.i, co2.i)
  # MNTD
  MNTD <- mntd.query(tree.i, co2.i, standardize = FALSE)
  MNTD.s <- mntd.query(tree.i, co2.i, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree.i, co2.i)

  out <- cbind(fips = row.names(co2.i), total.i, ontree.i, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("INTRO.county_diversity_INTROnull.", i-1, ".csv"), row.names = FALSE)
 
  # non-invasive introduced:
  # 1. null = all taxa in US regardless of status
  co2.ni <- co2
  co2.ni[, colnames(co2.ni)[!colnames(co2.ni) %in% ss.ni]] <- 0

  # PD
  PD <- pd.query(tree, co2.ni, standardize = FALSE)
  PD.s <- pd.query(tree, co2.ni, standardize = TRUE)
  PD.p <- pd.pvalues(tree, co2.ni)
  # MPD
  MPD <- mpd.query(tree, co2.ni, standardize = FALSE)
  MPD.s <- mpd.query(tree, co2.ni, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree, co2.ni)
  # MNTD
  MNTD <- mntd.query(tree, co2.ni, standardize = FALSE)
  MNTD.s <- mntd.query(tree, co2.ni, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree, co2.ni)

  out <- cbind(fips = row.names(co2.ni), total.ni, ontree.ni, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("NONINV.county_diversity_ALLnull.", i-1, ".csv"), row.names = FALSE)
  
  # 2. null = all NON-INVASIVE INTRODUCED taxa in US
  co2.ni <- co2[, colnames(co2) %in% ss.ni]  
  
  # PD
  PD <- pd.query(tree.ni, co2.ni, standardize = FALSE)
  PD.s <- pd.query(tree.ni, co2.ni, standardize = TRUE)
  PD.p <- pd.pvalues(tree.ni, co2.ni)
  # MPD
  MPD <- mpd.query(tree.ni, co2.ni, standardize = FALSE)
  MPD.s <- mpd.query(tree.ni, co2.ni, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree.ni, co2.ni)
  # MNTD
  MNTD <- mntd.query(tree.ni, co2.ni, standardize = FALSE)
  MNTD.s <- mntd.query(tree.ni, co2.ni, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree.ni, co2.ni)

  out <- cbind(fips = row.names(co2.ni), total.ni, ontree.ni, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("NONINV.county_diversity_NONINVnull.", i-1, ".csv"), row.names = FALSE)
  
  # 3. null = all INTRODUCED taxa in US regardless of invasive status
  co2.ni <- co2.i
  co2.ni[, colnames(co2.ni)[!colnames(co2.ni) %in% ss.ni]] <- 0
  
  # PD
  PD <- pd.query(tree.i, co2.ni, standardize = FALSE)
  PD.s <- pd.query(tree.i, co2.ni, standardize = TRUE)
  PD.p <- pd.pvalues(tree.i, co2.ni)
  # MPD
  MPD <- mpd.query(tree.i, co2.ni, standardize = FALSE)
  MPD.s <- mpd.query(tree.i, co2.ni, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree.i, co2.ni)
  # MNTD
  MNTD <- mntd.query(tree.i, co2.ni, standardize = FALSE)
  MNTD.s <- mntd.query(tree.i, co2.ni, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree.i, co2.ni)

  out <- cbind(fips = row.names(co2.ni), total.ni, ontree.ni, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("NONINV.county_diversity_INTROnull.", i-1, ".csv"), row.names = FALSE)
   
  # INVASIVE introduced:
  # 1. null = all taxa in US regardless of status
  co2.ii <- co2
  co2.ii[, colnames(co2.ii)[!colnames(co2.ii) %in% ss.ii]] <- 0

  # PD
  PD <- pd.query(tree, co2.ii, standardize = FALSE)
  PD.s <- pd.query(tree, co2.ii, standardize = TRUE)
  PD.p <- pd.pvalues(tree, co2.ii)
  # MPD
  MPD <- mpd.query(tree, co2.ii, standardize = FALSE)
  MPD.s <- mpd.query(tree, co2.ii, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree, co2.ii)
  # MNTD
  MNTD <- mntd.query(tree, co2.ii, standardize = FALSE)
  MNTD.s <- mntd.query(tree, co2.ii, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree, co2.ii)

  out <- cbind(fips = row.names(co2.ii), total.ii, ontree.ii, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("INV.county_diversity_ALLnull.", i-1, ".csv"), row.names = FALSE)
  
  # 2. null = all INVASIVE taxa in US
  co2.ii <- co2[, colnames(co2) %in% ss.ii]  
  
  # PD
  PD <- pd.query(tree.ii, co2.ii, standardize = FALSE)
  PD.s <- pd.query(tree.ii, co2.ii, standardize = TRUE)
  PD.p <- pd.pvalues(tree.ii, co2.ii)
  # MPD
  MPD <- mpd.query(tree.ii, co2.ii, standardize = FALSE)
  MPD.s <- mpd.query(tree.ii, co2.ii, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree.ii, co2.ii)
  # MNTD
  MNTD <- mntd.query(tree.ii, co2.ii, standardize = FALSE)
  MNTD.s <- mntd.query(tree.ii, co2.ii, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree.ii, co2.ii)

  out <- cbind(fips = row.names(co2.ii), total.ii, ontree.ii, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("INV.county_diversity_INVnull.", i-1, ".csv"), row.names = FALSE)
   
  # 3. null = all INTRODUCED taxa in US regardless of invasive status
  co2.ii <- co2.i
  co2.ii[, colnames(co2.ii)[!colnames(co2.ii) %in% ss.ii]] <- 0

  # PD
  PD <- pd.query(tree.i, co2.ii, standardize = FALSE)
  PD.s <- pd.query(tree.i, co2.ii, standardize = TRUE)
  PD.p <- pd.pvalues(tree.i, co2.ii)
  # MPD
  MPD <- mpd.query(tree.i, co2.ii, standardize = FALSE)
  MPD.s <- mpd.query(tree.i, co2.ii, standardize = TRUE)
  MPD.p <- mpd.pvalues(tree.i, co2.ii)
  # MNTD
  MNTD <- mntd.query(tree.i, co2.ii, standardize = FALSE)
  MNTD.s <- mntd.query(tree.i, co2.ii, standardize = TRUE)
  MNTD.p <- mntd.pvalues(tree.i, co2.ii)

  out <- cbind(fips = row.names(co2.ii), total.ii, ontree.ii, PD, PD.s, PD.p, MPD, MPD.s, MPD.p, MNTD, MNTD.s, MNTD.p)
  write.csv(out, paste0("INV.county_diversity_INTROnull.", i-1, ".csv"), row.names = FALSE)
  
}

stopCluster(cl)
gc()






