#Jamie R. Stavert #Centre for Biodiversity and Biosecurity #School of Biological Sciences #University of Auckland #Auckland, 1072 #New Zealand #12 June 2017 #jamie.stavert[at]gmail.com #This script assigns species to functional groups using weighted functional traits #It also returns deltaAICc values and evidence ratios for additive vs. interaction models for each functional group and additional species sets #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #load packages #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- require(cluster) require(NbClust) require(maptree) require(FD) require(plyr) require(dplyr) require(reshape2) require(ggplot2) require(MuMIn) require(glmmADMB) require(bbmle) require(ecodist) require(spdep) require(RANN) #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #read in data #"traits" contains values for functional traits #"response_divsersity" contains abundance for species at all sites #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- e <- read.csv("H:/Jamie/PhD Documents/Response diversity paper/ProcB Submission/Data & code/traits.csv", header=TRUE, row.names=1) a <-read.csv("H:/Jamie/PhD Documents/Response diversity paper/Response diversity spreadsheets/response_diversity_V2.csv", header=TRUE, row.names=NULL) #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #create trait weighting vector so that all trait grouping = 1 (hairiness; phenology; activity times; bodysize; tounge length; pollen collection structure) #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- w <- c(0.142857143, 0.142857143, 0.142857143, 0.142857143, 0.142857143, 0.142857143, 0.142857143,#hairiness 0.083333333, 0.083333333, 0.083333333, 0.083333333, 0.083333333, 0.083333333, 0.083333333, 0.083333333,0.083333333, 0.083333333, 0.083333333, 0.083333333,#phenology 0.142857143, 0.142857143, 0.142857143, 0.142857143, 0.142857143, 0.142857143, 0.142857143,#activity times 0.142857143, 0.142857143, 0.142857143, 0.142857143, 0.142857143, 0.142857143, 0.142857143,#body size 1,#tongue length 1)#pollen collection structure #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #create Gower dissimilarity matrix (with traits weighted) #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- e.dist <- gowdis(e, w) #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #find optimal number of functional groups #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- noclus <- agnes(e.dist, method="ward") b <- kgs(noclus, noclus$diss, maxclust=21)#6 clusters has lowest penalty score plot(names (b), b, xlab="Number of Clusters", ylab="Penalty score") #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #silhouette plot #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- pamx <- pam(e.dist, 6) summary(pamx) plot(pamx) #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #plot dendrogram #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- e.clust <- hclust(e.dist, method="ward.D2") plot(e.clust, main = "Cluster dengrogram based on effect traits") cut.g <- readline("6") cut.g <- as.integer(cut.g) e.gr <- cutree(e.clust, k = 6) e.gr2 <- rect.hclust(e.clust, k = 6, border = "red") #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #merge functional effect groups with species abundances dataframe #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- groups <- as.data.frame(e.gr) groups %>% tbl_df() groups <- add_rownames(groups, var = "species") a_groups <- merge(a, groups, by = "species", all.x = TRUE) a_groups$site <- as.factor(a_groups$site) #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for full community #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- options(na.action = "na.fail") summary(full_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data = a_groups, family="nbinom1"))#best model summary(full_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data = a_groups, family="nbinom2")) AICtab(full_int_a, full_int_b) plot(residuals(full_int_a,type=c("pearson"))) full_sel <- as.data.frame(dredge(full_int_a))#produces several warnings for abundance ~ (1 | site) and abundance ~ pagri + (1 | site) models #this is not a problem as we are not interested in these models full_sel full_sel[3,7]-full_sel[1,7]#calculate delta AICc full_sel[1,9]/full_sel[3,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for effect group 1 #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- summary(eg1_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="1"), family="nbinom1"))#best model summary(eg1_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="1"), family="nbinom2")) AICtab(eg1_int_a, eg1_int_b) plot(residuals(eg1_int_a,type=c("pearson"))) eg1_sel <- as.data.frame(dredge(eg1_int_a)) eg1_sel eg1_sel[1,7]-eg1_sel[3,7]#calculate delta AIC eg1_sel[3,9]/eg1_sel[1,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for effect group 2 #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- summary(eg2_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="2"), family="nbinom1")) summary(eg2_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="2"), family="nbinom2"))#best model AICtab(eg2_int_a, eg2_int_b) plot(residuals(eg2_int_b,type=c("pearson"))) eg2_sel <- as.data.frame(dredge(eg2_int_b)) eg2_sel eg2_sel[2,7]-eg2_sel[3,7]#calculate delta AIC eg2_sel[3,9]/eg2_sel[2,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for effect group 3 #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- summary(eg3_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="3"), family="nbinom1"))#unable to estimate standard errors with this summary(eg3_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="3"), family="nbinom2"))#best model plot(residuals(eg3_int_b,type=c("pearson"))) eg3_sel <- as.data.frame(dredge(eg3_int_b)) eg3_sel eg3_sel[1,7]-eg3_sel[3,7]#calculate delta AIC eg3_sel[3,9]/eg3_sel[1,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for effect group 4 #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- summary(eg4_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="4"), family="nbinom1"))#best model summary(eg4_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="4"), family="nbinom2")) AICtab(eg4_int_a, eg4_int_b) plot(residuals(eg4_int_a,type=c("pearson"))) eg4_sel <- as.data.frame(dredge(eg4_int_a)) eg4_sel eg4_sel[2,7]-eg4_sel[1,7]#calculate delta AIC eg4_sel[1,9]/eg4_sel[2,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for effect group 6 #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- summary(eg6_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="6"), family="nbinom1"))#best model summary(eg6_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,e.gr=="6"), family="nbinom2")) AICtab(eg6_int_a, eg6_int_b) plot(residuals(eg6_int_a,type=c("pearson"))) eg6_sel <- as.data.frame(dredge(eg6_int_a))#unable to estimate standard errors with this for abundance ~ (1 | site) model #not a problem as we don't require this model eg6_sel eg6_sel[4,7]-eg6_sel[1,7]#calculate delta AIC eg6_sel[1,9]/eg6_sel[4,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for exotics #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- summary(exotic_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,exotic_native=="exotic"), family="nbinom1"))#best model summary(exotic_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,exotic_native=="exotic"), family="nbinom2")) AICtab(exotic_int_a, exotic_int_a) plot(residuals(exotic_int_a,type=c("pearson"))) exotic_sel <- as.data.frame(dredge(exotic_int_b))#convergence failed for abundance ~ (1 | site) model #not a problem as we don't require this model exotic_sel exotic_sel[2,7]-exotic_sel[1,7]#calculate delta AIC exotic_sel[1,9]/exotic_sel[2,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for natives #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- summary(native_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,exotic_native=="native"), family="nbinom1")) summary(native_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,exotic_native=="native"), family="nbinom2"))#best model AICtab(native_int_a, native_int_b) plot(residuals(native_int_b,type=c("pearson"))) native_sel <- as.data.frame(dredge(native_int_b)) native_sel native_sel[2,7]-native_sel[1,7]#calculate delta AIC native_sel[1,9]/native_sel[2,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for bees #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- summary(bee_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,bee_fly=="bee"), family="nbinom1"))#best model summary(bee_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,bee_fly=="bee"), family="nbinom2")) AICtab(bee_int_a, bee_int_b) plot(residuals(bee_int_a,type=c("pearson"))) bee_sel <- as.data.frame(dredge(bee_int_a)) bee_sel bee_sel[3,7]-bee_sel[1,7]#calculate delta AIC bee_sel[1,9]/bee_sel[3,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #run model and calculate delta AICc for flies #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- summary(fly_int_a <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,bee_fly=="fly"), family="nbinom1"))#best model summary(fly_int_b <- glmmadmb(abundance ~ pagri * species + (1|site), data=subset(a_groups,bee_fly=="fly"), family="nbinom2")) AICtab(fly_int_a, fly_int_b) plot(residuals(fly_int_a,type=c("pearson"))) fly_sel <- as.data.frame(dredge(fly_int_a)) fly_sel fly_sel[3,7]-fly_sel[1,7]#calculate delta AIC fly_sel[1,9]/fly_sel[3,9]#calculate evidence ratio #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #check for spatial autocorrelation #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #read in coordinates xy <- read.csv("H:\\Jamie\\PhD Documents\\Response diversity paper\\Response diversity spreadsheets\\Landuse\\site_coordinates.csv",header=TRUE, row.names=NULL) #run moran test for spatial autocorrelation d <- merge(a, xy[,-2], by = "site", all.x = TRUE) xy <- dplyr::select(d, one_of(c("lat", "long"))) xy <- as.matrix(xy) UTMnb <- knn2nb(knearneigh(jitter(xy), k=nrow(xy)-1)) UTMls <- nb2listw(UTMnb) #make a list moran(d$abundance, UTMls, length(UTMnb), Szero(UTMls)) #Moran index for data rd <- glmmadmb(abundance ~ pagri * species + (1|site), data = a_groups, family="nbinom1") moran(residuals(rd, type="response"), UTMls, length(UTMnb), Szero(UTMls))#Moran for residuals ############################ #END ############################