###Text S1.R Code produce most statistics in this publication broken into 4 parts: 1) Null simulations using DAMOCLES, 2) Phylogenetic community turnover, 3) phylogenetic community outlier detection, and 4) GLMM models predicting island presence given mainland abundance. #### See Harris et al. 2015 for NGMS package for UNTB simulations ####Make sure you have all of the datasets in a single file for this adn set your directory to that file. library(picante) library(DAMOCLES) library(ape) ###PART 1 ###This first set of codes are essentially a wrapper for Pigot and Etienne's (2015) DAMOCLES. It will take your communities, in my case I did one set for islands and one set ###for the adjacent mainland, estimate colonization (gamma) and mu (extinction), simulate null models on those and compare your real values estimated from various community ###metrics such as psv, mntd, and mpd to the simulated values. ###Make sure your phylogenies and communities are matched prior to using this (I used match.phy.comm in picante) ###Run these example data from my paper. read.table("ISL_PRES", header=T, row.names=1)->isl read.tree("tree.txt")->tree match.phylo.comm(phy=tree, comm=t(isl))->matched2 ###The following will generate damocles parameters gamma and mu for your communities, see below for running ###everything. Note, you can change the defaults in damocles easily ###by providing different parameters for these functions goddamocles<-function(tree, comm, initparsopt, idparsopt, parsfix, idparsfix, pars2, pchoice) {lapply(c(1:length(row.names(comm))), function(x) cbind(row.names(t(comm)), as.matrix(as.numeric(comm[x,]))))->comm2; lapply(comm2, function (x) DAMOCLES_ML(phy = tree,pa = x, initparsopt, idparsopt, parsfix, idparsfix,pars2,pchoice))->out; return(out) } ###Run it the following way on a matched object (community matched with phylogeny) from picante (named matched2 here). IMPORTANT: I made an island community matrix that ##has all US species even if they don't appear on the islands ISL_PRES and matched this to the tree of 144 taxa. I did this so that I could use all of the taxa ###from the US or speciation info and put those taxa not appearing in any of the islands into the island communit matrix with "0" listed for all of the localities of course. I then ##matched this new island comm to the phylogeny of all taxa. This was again performed for the mainland communities. All mainland and parallel island ###communities are sorted ##the same in the two files. goddamocles(matched2$phy, comm=matched2$comm, initparsopt=c(0.05,0.05), idparsopt=c(1,2), parsfix=c(0), idparsfix=c(3), pars2=c(1E-3,1E-4,1E-5,1000), pchoice = 0)->output ###get all of your desired gamma and mu outputs for each island from the output via lapply(c(1:length(output)), function(x) output[[x]]$gamma_0)->gamma lapply(c(1:length(output)), function(x) output[[x]]$mu)->mu ###Here are the functions to perform the simulations from Pigot & Etienne (2015) from the first community using their estimated respective mu and gamma, again you can change the initial params however. DAMOCLES_sim2<- function (phy, gamma_0, gamma_td, mu, sigma, psiBranch, psiTrait, z, phi, traitOpt, br0, br_td, nTdim, root.state, root.trait.state, plotit = FALSE, keepExtinct = FALSE) { if (plotit == TRUE) { plot(phy) } dn = round(dist.nodes(phy), 4) ntips = length(phy$tip.label) nbranch = 2 * ntips - 2 patable = data.frame(p = phy$edge[, 1], d = phy$edge[, 2], age.start = rep(NA, nbranch), age.end = rep(NA, nbranch), extant = rep(0, nbranch), state = rep(NA, nbranch)) patable$age.start = dn[patable$p, ntips + 1] patable$age.end = dn[patable$d, ntips + 1] ce = which(patable$p == min(patable$p)) patable$extant[ce] = 1 patable$state = rep(NA, dim(patable)[1]) patable$state[ce] = sample(c(0, root.state)) if (nTdim > 1) { traits = matrix(NA, ncol = dim(patable)[1], nrow = nTdim) traits[, ce] = 0 } else { traits = rep(NA, dim(patable)[1]) traits[ce] = 0 } tstep = min(patable$age.start) while (tstep < max(patable$age.end)) { extant = which(patable$extant == 1) extant.id = patable$d[extant] pa = patable$state[extant] Itot.max = gamma_0 * length(pa) Etot.max = mu * length(pa) tstep0 = tstep wt = -log(runif(1))/(Itot.max + Etot.max) tstep1 = tstep0 + wt tstep = min(tstep1, min(patable$age.end[extant])) if (br_td == 0) { wt = tstep - tstep0 } else { wt = (max(patable$age.end) * (1 - exp(br_td * tstep))) - (max(patable$age.end) * (1 - exp(br_td * tstep0))) } sdnorm = sqrt(br0 * wt) traits = traits + rnorm(length(traits), mean = 0, sd = sdnorm) if (tstep < max(patable$age.end)) { if (tstep < min(patable$age.end[extant])) { Itot.real = 0 if (length(pa) - sum(pa) > 0) { Dpaynotpa = 1 - pa Itot.real = gamma_0 * sum(Dpaynotpa) if (gamma_td > 0) { Dpaynotpa = 1 - pa Dpaynotpa = gamma_0/(1 + (gamma_td * tstep)) * Dpaynotpa Itot.real = sum(Dpaynotpa) } if (psiBranch > 0 & sum(pa) > 0) { Db = matrix(rep(patable$age.end[extant] - tstep, sum(patable$extant)), nrow = sum(patable$extant)) futD = Db + t(Db) D = dn[extant.id, extant.id] - futD D = replace(D, D < 0, 0) D[, which(pa == 1)] = NaN D[which(pa == 0), ] = NaN Dpay = D^z/(psiBranch^z + D^z) Dpaynotpa = colProds(Dpay, na.rm = TRUE) * (1 - pa) Dpaynotpa = gamma_0 * Dpaynotpa Itot.real = sum(Dpaynotpa) } if (phi > 0) { if (nTdim > 1) { D = traitOpt - traits[, extant] D = sqrt(colSums(D^2)) } else { D = abs(traitOpt - traits[extant]) } Dpaynotpa = dexp(D, rate = phi)/dexp(0, rate = phi) Dpaynotpa = gamma_0 * Dpaynotpa Dpaynotpa = Dpaynotpa * (1 - pa) Itot.real = sum(Dpaynotpa) } if (psiTrait > 0 & sum(pa) > 0) { if (nTdim > 1) { D <- dist(t(traits[, extant]), method = "euclidean", diag = FALSE, upper = TRUE, p = 2) } else { D <- dist(traits[extant], method = "euclidean", diag = FALSE, upper = TRUE, p = 2) } D = as.matrix(D) D[which(pa == 0), ] = NaN Dpay = D^z/(psiTrait^z + D^z) Dpaynotpa = colProds(Dpay, na.rm = TRUE) * (1 - pa) Dpaynotpa = gamma_0 * Dpaynotpa Itot.real = sum(Dpaynotpa) } } Etot.real = mu * sum(pa) Itot = Itot.real/(Itot.max + Etot.max) Etot = Etot.real/(Itot.max + Etot.max) event = sample(x = 0:2, size = 1, replace = F, prob = c(1 - (Itot + Etot), Itot, Etot)) if (event == 1) { if (sum(1 - pa) == 1) { patable$state[extant[which(pa == 0)]] = 1 } else { patable$state[sample(extant, 1, replace = FALSE, prob = Dpaynotpa)] = 1 } if (plotit == TRUE) { lines(c(tstep, tstep), c(0, 100), col = "lightblue", lty = 2) } } if (event == 2) { if (sum(pa) == 1) { patable$state[extant[which(pa == 1)]] = 0 } else { patable$state[sample(extant[which(pa == 1)], 1)] = 0 } if (plotit == TRUE) { lines(c(tstep, tstep), c(0, 100), col = "orange", lty = 2) } } } else { focedge = extant[which(patable$age.end[extant] == min(patable$age.end[extant]))] for (fe in unique(focedge)) { dedge = which(patable$p == patable$d[fe]) if (length(dedge) > 0) { patable$extant[fe] = 0 patable$extant[dedge] = 1 if (nTdim > 1) { traits[, dedge] = traits[, fe] } else { traits[dedge] = traits[fe] } ds = sample(c(0, patable$state[fe]), 1, prob = c(1 - sigma, sigma)) patable$state[dedge] = sample(c(ds, patable$state[fe])) } else { patable$extant[fe] = 0 } } } } } patablelist = list() if (keepExtinct == FALSE) { if (nTdim > 1) { traits = traits[, which(patable$extant == 1)] } else { traits = traits[which(patable$extant == 1)] } patable = patable[which(patable$extant == 1), ] } patablelist[[1]] = patable #patablelist[[2]] = traits return(patablelist) } ###Now, run this to get your simulated null communities given gamma and mu, added below in the function as output[[x]]$gamma_0 and output[[x]]$mu. This can take a while. DAM_ISL_SIM <- lapply(c(1:length(output)), function(x) replicate(100, (DAMOCLES_sim2( phy=matched2$phy, gamma_0 = output[[x]]$gamma_0, gamma_td =0, mu = output[[x]]$mu, sigma = 0, psiBranch = 0, psiTrait = 0, z = 0, phi = 0, traitOpt = 0, br0 = 0, br_td = 0, nTdim = 0, root.state = 1, root.trait.state = 0, plotit = FALSE, keepExtinct = FALSE ))[[1]]$state, FALSE)) ###this will get you the P values for overdispersion for the community sims from damocles vs. your real community, assume you made 100 replicates using the phylogenetic community metrics psv, or mntd, mpd. Here comm_psv<-function(sim, tree, comm) {lapply(sim, function(x) as.data.frame(x,row.names=tree$tip.label))->nulls lapply(nulls, function(x) match.phylo.comm(tree, t(x)))->matched lapply(c(1:length(matched)), function(x) psv(matched[[x]]$comm, matched[[x]]$phy)[[1]])->simPSVs match.phylo.comm(tree, t(comm))->matchreal psv(matchreal$comm, matchreal$phy)$PSVs->realPSV lapply(c(1:length(simPSVs)), function(x) mean(simPSVs[[x]], na.rm=T))->sim_means lapply(c(1:length(simPSVs)), function(x) sd(simPSVs[[x]], na.rm=T))->sim_sd unlist(lapply(c(1:length(sim_means)), function(x) (realPSV[x]-sim_means[[x]])/sim_sd[[x]]))->zscores (unlist(lapply(c(1:length(simPSVs)), function(x) rank(c(realPSV[x], simPSVs[[x]]))[1])))->PSV.rank PSV.obs.rank <- ifelse(is.na(zscores), NA, PSV.rank) PSV.obs.rank /(length(simPSVs[[1]])+1)->p_val results<-list() results$Z_scores<-zscores results$P<-p_val return(results) } comm_mntd<-function(sim, tree, comm) {lapply(sim, function(x) as.data.frame(x,row.names=tree$tip.label))->nulls lapply(nulls, function(x) match.phylo.comm(tree, t(x)))->matched lapply(c(1:length(matched)), function(x) mntd(matched[[x]]$comm, cophenetic(matched[[x]]$phy)))->simmntd match.phylo.comm(tree, t(comm))->matchreal mntd(matchreal$comm, cophenetic(matchreal$phy))->realmntd lapply(c(1:length(simmntd)), function(x) mean(simmntd[[x]], na.rm=T))->sim_means lapply(c(1:length(simmntd)), function(x) sd(simmntd[[x]], na.rm=T))->sim_sd unlist(lapply(c(1:length(sim_means)), function(x) (realmntd[x]-sim_means[[x]])/sim_sd[[x]]))->zscores (unlist(lapply(c(1:length(simmntd)), function(x)rank(c(realmntd[x], simmntd[[x]]))[1])))->mntd.rank mntd.obs.rank <- ifelse(is.na(zscores), NA, mntd.rank) mntd.obs.rank /(length(simmntd[[1]])+1)->p_val results<-list() results$Z_scores<-zscores results$P<-p_val return(results) } comm_mpd<-function(sim, tree, comm) {lapply(sim, function(x) as.data.frame(x,row.names=tree$tip.label))->nulls lapply(nulls, function(x) match.phylo.comm(tree, t(x)))->matched lapply(c(1:length(matched)), function(x) mpd(matched[[x]]$comm, cophenetic(matched[[x]]$phy)))->simmpd match.phylo.comm(tree, t(comm))->matchreal mpd(matchreal$comm, cophenetic(matchreal$phy))->realmpd lapply(c(1:length(simmpd)), function(x) mean(simmpd[[x]], na.rm=T))->sim_means lapply(c(1:length(simmpd)), function(x) sd(simmpd[[x]], na.rm=T))->sim_sd unlist(lapply(c(1:length(sim_means)), function(x) (realmpd[x]-sim_means[[x]])/sim_sd[[x]]))->zscores (unlist(lapply(c(1:length(simmpd)), function(x)rank(c(realmpd[x], simmpd[[x]]))[1])))->mpd.rank mpd.obs.rank <- ifelse(is.na(zscores), NA, mpd.rank) mpd.obs.rank /(length(simmpd[[1]])+1)->p_val results<-list() results$Z_scores<-zscores results$P<-p_val return(results)} ###For example, using the mpd stat we would do the following: ###So for islands, remember from what we have done far for these functions sim=damocles sim, ####tree=phylogenetic tree, comm=isl. read.table("Appendix I.Sites_Presence_Abundance.txt", header=T, row.names=1)->sites as.matrix(sites[seq(1,162,3)])->isl as.matrix(sites[seq(2,162,3)])->main ###you can simulate the island, or pull the file here: dget("DAM_ISL_SIM")->DAM_ISL_SIM comm_psv(sim= DAM_ISL_SIM,tree= tree,comm= isl)->outPSV_isl comm_mntd(sim= DAM_ISL_SIM,tree= tree,comm= isl)->outmntd_isl comm_mpd(sim= DAM_ISL_SIM,tree= tree,comm= isl)->outMPD_isl ###you can simulate the mainland as above with the island, or pull the file here: dget("DAM_MAIN_SIM")->DAM_MAIN_SIM comm_psv(sim= DAM_MAIN_SIM,tree= tree,comm= main)->outPSV_main comm_mntd(sim= DAM_MAIN_SIM,tree= tree,comm= main)->outmntd_main comm_mpd(sim= DAM_MAIN_SIM,tree= tree,comm= main)->outMPD_main ###END of PART1 ###PART 2, this will compare intercommunity phylpogenetic distance using a few different metrics. I use the same tree as above but not a new file that includes the island ###and mainlaind community presence and abundance. We will only use presence. First we need to procure our island and mainland files. Mine are sewn together ###into a single ###file. read.table("Appendix I.Sites_Presence_Abundance.txt", header=T, row.names=1)->sites as.matrix(sites[seq(1,162,3)])->ISL as.matrix(sites[seq(2,162,3)])->MAIN ###This gets you the phylogenetic distance between each paired island-mainland community, it takes a bit of time to finish. PCDpshuffle<-function(main, isl, tree,reps) { as.matrix(main)->main as.matrix(isl)->isl sum(main)-sum(isl)->diff which(main==1)->list replicate(reps,sample(list, diff, T))->tests which(main==1)->start as.matrix(start)->start1 start1[,1]<-1 cbind(start1,replicate(reps,rep(1,length(start1))))->newStart replacer<-function(comm, diff) { comm[sample(1:length(comm), diff)]<-0 return(comm) } lapply(c(1:reps), function(x) replacer(newStart[,x], diff))->out lapply(c(1:reps), function(x) as.numeric(out[[x]]))->out1 as.data.frame(out1)->outw cbind(start1, outw)->test colnames(test)<-c(1:(reps+1)) row.names(test)<-row.names(main)[list] match.phylo.comm(tree, t(test))->match pcd(t(test),tree, reps=reps)$PCDp->junk as.matrix(junk)[,1][-1]->dist mean(as.matrix(junk)[,1][-1], na.rm=T)->meanDist sd(as.matrix(junk)[,1][-1], na.rm=T)->sdDist cbind (main, isl)->real real[ rowSums(real==0) != ncol(real), ] ->real pcd(t(real), tree, reps=reps)$PCDp->stat (stat-meanDist)/sdDist->z_score rank(c(stat, dist))[1]->rankV rankV/(length(dist)+1)->p_val results<-list() results$stat<-stat results$Z_score<-z_score results$p_val<-p_val return(results) } ###Run this over all of your datasets (remember to make you communities matrices). PCD takes a long time to do 1000 reps, so I just do 10 of 100, but you could set this for 1000 straightaway. Also, I have 54 paired communities, you should change that for the number of communities you acually have. replicate(10,lapply(c(1:54), function(x) PCDpshuffle(MAIN[,x],ISL[,x], tree=tree, reps=100)))->outPCD ### I ran this using the replicate function becuase it takes forever to do 1000 runs, so replicate at 10. But the output is indexed horribly, so use: get_out<-function(dataset, number, index, type) { sapply(dataset[number, ][index], "[[", type)->out return(out) } ###Get the P and Z score here lapply(c(1:54), function(x) mean(get_out(dataset=outPCD, number=x, index=c(1:10), "p_val")) )->PCD_P lapply(c(1:54), function(x) mean(get_out(dataset=outPCD, number=x, index=c(1:10), "Z_score")) )->PCD_Z ###Does the same as above but for mean nearest taxon distance mntdshuffle<-function(main, isl, tree, reps) { as.matrix(main)->main as.matrix(isl)->isl sum(main)-sum(isl)->diff which(main==1)->list replicate(reps,sample(list, diff, T))->tests which(main==1)->start as.matrix(start)->start1 start1[,1]<-1 cbind(start1,replicate(reps,rep(1,length(start1))))->newStart replacer<-function(comm, diff) { comm[sample(1:length(comm), diff)]<-0 return(comm) } lapply(c(1:reps), function(x) replacer(newStart[,x], diff))->out lapply(c(1:reps), function(x) as.numeric(out[[x]]))->out1 as.data.frame(out1)->outw cbind(start1, outw)->test colnames(test)<-c(1:(reps+1)) row.names(test)<-row.names(main)[list] match.phylo.comm(tree, t(test))->match comdistnt(comm=match$comm, dis=cophenetic(match$phy), abundance.weighted = FALSE, exclude.conspecifics = FALSE)[1:reps]->dist cbind (main, isl)->real real[ rowSums(real==0) != ncol(real), ] ->real match.phylo.comm(tree, t(real))->matchR comdistnt(comm=matchR$comm, dis=cophenetic(matchR$phy), abundance.weighted = FALSE, exclude.conspecifics = FALSE)[1]->stat mean(dist, na.rm=T)->meanDist sd(dist, na.rm=T)->sdDist (stat-meanDist)/sdDist->z_score rank(c(stat, dist))[1]->rankV rankV/(length(dist)+1)->p_val results<-list() results$Z_score<-z_score results$p_val<-p_val return(results)} lapply(c(1:54), function(x) mntdshuffle(MAIN[,x],ISL[,x], tree=tree, reps=100))->outMNTDshuffle ###get the prob and Z score sapply(c(1:54), function(x) outMNTDshuffle[[x]]$p_val)->outMNTDshuffleP sapply(c(1:54), function(x) outMNTDshuffle[[x]]$Z_score)->outMNTDshuffleZ ###does the same as above but for mpd MPDshuffle<-function(main, isl, tree, reps) { as.matrix(main)->main as.matrix(isl)->isl sum(main)-sum(isl)->diff which(main==1)->list replicate(reps,sample(list, diff, T))->tests which(main==1)->start as.matrix(start)->start1 start1[,1]<-1 cbind(start1,replicate(reps,rep(1,length(start1))))->newStart replacer<-function(comm, diff) { comm[sample(1:length(comm), diff)]<-0 return(comm) } lapply(c(1:reps), function(x) replacer(newStart[,x], diff))->out lapply(c(1:reps), function(x) as.numeric(out[[x]]))->out1 as.data.frame(out1)->outw cbind(start1, outw)->test colnames(test)<-c(1:(reps+1)) row.names(test)<-row.names(main)[list] match.phylo.comm(tree, t(test))->match 1-comdist(comm=match$comm, dis=cophenetic(match$phy), abundance.weighted = FALSE)[1:reps]->dist cbind (main, isl)->real real[ rowSums(real==0) != ncol(real), ] ->real match.phylo.comm(tree, t(real))->matchR 1-comdist(comm=matchR$comm, dis=cophenetic(matchR$phy), abundance.weighted = FALSE)[1]->stat mean(dist, na.rm=T)->meanDist sd(dist, na.rm=T)->sdDist (stat-meanDist)/sdDist->z_score rank(c(stat, dist))[1]->rankV rankV/(length(dist)+1)->p_val results<-list() results$Z_score<-z_score results$p_val<-p_val return(results)} ###run it lapply(c(1:54), function(x) MPDshuffle(MAIN[,x],ISL[,x], tree=tree, reps=100))->outMPDshuffle ###get the P and Z score sapply(c(1:54), function(x) outMPDshuffle[[x]]$p_val)->outMPDshuffleP sapply(c(1:54), function(x) outMPDshuffle[[x]]$Z_score)->outMPDshuffleZ ###Does the same as above but uses the metric from the Ricotta et al. 2015 called Dab DabShuffle<-function(main, isl, tree, reps) { as.matrix(main)->main as.matrix(isl)->isl sum(main)-sum(isl)->diff which(main==1)->list replicate(reps,sample(list, diff, T))->tests which(main==1)->start as.matrix(start)->start1 start1[,1]<-1 cbind(start1,replicate(reps,rep(1,length(start1))))->newStart replacer<-function(comm, diff) { comm[sample(1:length(comm), diff)]<-0 return(comm) } lapply(c(1:reps), function(x) replacer(newStart[,x], diff))->out lapply(c(1:reps), function(x) as.numeric(out[[x]]))->out1 as.data.frame(out1)->outw cbind(start1, outw)->test colnames(test)<-c(1:(reps+1)) row.names(test)<-row.names(main)[list] match.phylo.comm(tree, t(test))->match phylo_dist<-function(abundances, dissimilarities) { if(any(dissimilarities>1)){ warning("Phylogenetic dissimilarities are not in the range 0-1. They have been normalized by the maximum") dissimilarities <- dissimilarities/max(dissimilarities) } if(any(!colnames(abundances) %in%rownames(dissimilarities))) stop("At least one species in the matrix of abundances is missing in the matrix of dissimilarities") if(any(!colnames(abundances) %in%colnames(dissimilarities))) stop("At least one species in the matrix of abundances is missing in the matrix of dissimilarities") dissimilarities <- dissimilarities[colnames(abundances), colnames(abundances)] dataset <- t(abundances) similarities<-1-as.matrix(dissimilarities) total <- colSums(dataset) rel_abu<- sweep(dataset, 2, total, "/") num.plot<-dim(dataset)[2] names<-list(colnames(dataset), colnames(dataset)) dist.matrix<-matrix(0, nrow=num.plot, ncol=num.plot, dimnames=names) for (i in 2:num.plot) { for (j in 1:(i-1)) { mat_folk<-similarities*rel_abu[,i] mat_folk2<- similarities*rel_abu[,j] xxx<-colSums(mat_folk2) yyy<-colSums(mat_folk) sub<-abs(xxx-yyy) index<-sum(sub)/sum(xxx+yyy) dist.matrix[i,j]<-index } } dist.matrix <- dist.matrix + t(dist.matrix) semi_matrix<-as.dist(dist.matrix) return(semi_matrix) } as.matrix(phylo_dist(abundances=match$comm, dissimilarities=cophenetic(match$phy)))[,1][-1]->dist cbind (main, isl)->real real[ rowSums(real==0) != ncol(real), ] ->real match.phylo.comm(tree, t(real))->matchR phylo_dist(abundances=matchR$comm,dissimilarities=cophenetic(matchR$phy))[[1]]->stat mean(dist, na.rm=T)->meanDist sd(dist, na.rm=T)->sdDist (stat-meanDist)/sdDist->z_score rank(c(stat, dist))[1]->rankV rankV/(length(dist)+1)->p_val results<-list() results$Z_score<-z_score results$p_val<-p_val return(results)} lapply(c(1:54), function(x) DabShuffle(MAIN[,x],ISL[,x], tree=tree, reps=100))->outDab ###get the P and Z score sapply(c(1:54), function(x) outDab[[x]]$p_val)->outDabshuffleP sapply(c(1:54), function(x) outDab[[x]]$Z_score)->outDabshuffleZ ####PART3 ###determines if you have outliers that heavily influence phylogenetic diversity when comparing to localities (mainland and island), it will not work if you only have one individual in a community, so skip those or use my loop below library(outliers) library(ape) library(picante) library(extremevalues) PD_Impact<-function(dat_isl,dat_main,dat_tot, tree) { sum(dat_main)-sum(dat_isl)->diff which(dat_main==1)->list replicate(100,sample(list, diff, T))->tests #replicate(84,0)->val replicate(100, replicate(84,0))->val replacer<-function(val, tests) {val[tests]<-1 names(val)<-row.names(dat) return(val)} ###determine which individuals impact PSV the most, start with matched communities dat_isl->junky names(junky)<-row.names(dat_tot) as.matrix(junky)->junky cbind(junky, junky[,1])->trial match.phylo.comm(phy=tree, comm=t(trial))->match drop.tip (match$phy, which(match$comm[1,]==0))->dropped_tree lapply(c(1:length(dropped_tree$tip.label)), function(x) drop.tip(dropped_tree, dropped_tree$tip.label[x]))->testers make_comm<-function(tree) {replicate(length(tree$tip.label),1)->this names(this)<-tree$tip.label return(this)} lapply(testers, make_comm)->comms mapply(function (x, y) psv(x, y)[1], x=comms, y=testers)->out make_comm(dropped_tree)->dcomm psv(dcomm, dropped_tree)[1]->actual out-actual->stuffy stuffy/sum(abs(stuffy))->junk names(junk)<-dropped_tree$tip.label #sort(abs(junk))->junk2 sort(junk)->junk2 #grubbs.test(junk2, type=20, opposite=T)->try getOutliers(junk2, distribution="normal", method="II", alpha=c(0.025, 0.025))$iRight->try results<-list() junk2[try]->this results$value<-junk2 results$outliers<-this return(results)} ###Use this to loop it over a large community structure, some with only 1 taxon that would give errors naturally, and we catch them here. Here the structure is in the seq part such that you are comparing 1 (island) to 2 (mainland) and 4 (island) to mainlaind (5) etc. and catching errors that ##exist #for those with only 1 or two taxa. ###Here we are using sites from the my file, which the columns are Island_pres, Mainland_pres, and Mainland_abundance replicate 54 times. We only need the first two: read.tree("tree.txt")->tree read.table("Appendix I.Sites_Presence_Abundance.txt", header=T, row.names=1)->sites mapply(function(x,y) tryCatch({PD_Impact(sites[,x],sites[,y],sites, tree)},error=function(e)(NA)),seq(1,162,3), seq(2,162,3)) ####PART 4 ####Using the GLMM to predict presence or absence on islands library(lme4) #Get the datafile which is stored in R format and lists mainland abundance, presence on island and name of each insland along three columns dget("Island_Mainland_Freq" )->freq glmer(presence_island~log(freq_mainland)+(1|island)+ (1|names), data=freq, family=binomial ("logit"),control = glmerControl(optimizer = "Nelder_Mead"), nAGQ = 1)->glmOUT ###Get confidence intervals for the fixed effects from the SE se <- sqrt(diag(vcov(glmOUT))) (outCI<- cbind(Est = fixef(glmOUT), LL = fixef(glmOUT) - 1.96 * se, UL = fixef(glmOUT) + 1.96 * se)) ###as odds ratios exp(outCI)