---
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")
```
|