--- title: "R Notebook" output: html_notebook --- ```{r} rm(list = ls()) library(indicspecies) library(data.table) set.seed(250) ``` #Read in invertebrate community data and read in and organize invertebrate sample metadata ```{r} df <- read.csv("~/Dropbox/Active/Arizona/1) CA Smith Fellowship/Ecosphere Submission/Datasets/DataS3.csv", header = T, row.names = 1, skip=1) df <- df[,c(7:88)] df <- t(df) md <- read.csv("~/Dropbox/Active/Arizona/1) CA Smith Fellowship/Ecosphere Submission/Datasets/DataS2.csv", header = T, skip=1) md <- md[md$Sampling.x.site.code!="",] row.names(md) <- md$Sampling.x.site.code md <- md[,-1] md <- md[rownames(df), ] md$Season <- as.factor(sapply(strsplit(as.character(md$Season.year)," "), `[`, 1)) ``` ###Overall test of indicator species (by basin and flow class) ```{r} md$Basin.Flow <- paste(md$Basin, md$Flow.class, sep=" ") indvalrest = multipatt(df, md$Basin.Flow, func = "r.g", control = how(nperm=1000)) indvalrest.sign<-as.data.table(indvalrest$sign, keep.rownames=TRUE) indvalrest.sign[ ,p.value.bh:=p.adjust(p.value, method="BH")] res <- indvalrest.sign[p.value.bh<=0.05, ] write.csv(res, "AllDataInvertIndicator.csv") ``` ###Subset by basin ```{r} md.pin <- md[which(md$Basin=="Pinnacles"),] df.pin <- df[c(row.names(md.pin)),] row.names(md.pin) == row.names(df.pin) md.pr<- md[which(md$Basin=="Point Reyes"),] df.pr <- df[c(row.names(md.pr)),] row.names(md.pr) == row.names(df.pr) ``` ###Indicators of flow class and season Pinnacles ```{r} md.pin$Flow.Season <- paste(md.pin$Flow.class, md.pin$Season, sep=" ") indvalrest = multipatt(df.pin, md.pin$Flow.Season, func = "r.g", control = how(nperm=1000)) indvalrest.sign<-as.data.table(indvalrest$sign, keep.rownames=TRUE) indvalrest.sign[ ,p.value.bh:=p.adjust(p.value, method="BH")] res <- indvalrest.sign[p.value.bh<=0.05, ] write.csv(res, "PinInvertIndicator.csv") ``` ###Indicators of flow class and season Point Reyes ```{r} md.pr$Flow.Season <- paste(md.pr$Flow.class, md.pr$Season, sep=" ") indvalrest = multipatt(df.pr, md.pr$Flow.Season, func = "r.g", control = how(nperm=1000)) indvalrest.sign<-as.data.table(indvalrest$sign, keep.rownames=TRUE) indvalrest.sign[ ,p.value.bh:=p.adjust(p.value, method="BH")] res <- indvalrest.sign[p.value.bh<=0.05, ] write.csv(res, "PrInvertIndicator.csv") ``` |