model_native_nonnative.Rmd

---
title: "model_native_nonnative"
author: "NCEAS"
date: '2022-08-11'
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
library(readr)
library(ggplot2)
library(dplyr)
library(tidyr)

library(lme4)
library(lmerTest)
library(emmeans)

library(sjPlot)
```

# Read data
```{r}
data_exceedance0 <- read_csv("SpeciesTables/daysExceedanceUpperDataForPlotting.csv") %>%
  filter(!Species %in% c("Giant Reed EAV", "Brazilian Waterweed", "Water Hyacinth FAV", "Microcystis", "Asian Clam", "Overbite Clam"))
```

# Categorize data
```{r}
native <- c("Green Sturgeon", "White Sturgeon","Delta Smelt", "Longfin Smelt", "Wakasagi","Chinook Salmon", "Steelhead/Rainbow Trout","Sacramento Splittail", "Sacramento Pikeminnow", "Sacramento Sucker", "Sacramento Blackfish", "Hitch", "Tule Perch", "Threespine Stickleback", "Prickly Sculpin")

# nonnative <- c("Mississippi Silverside", "Striped Bass", "Largemouth Bass", "Giant Reed EAV", "Brazilian Waterweed", "Water Hyacinth FAV", "Microcystis", "Asian Clam", "Overbite Clam")
nonnative <- c("Mississippi Silverside", "Striped Bass", "Largemouth Bass")

listed <- c("Green Sturgeon","Delta Smelt", "Longfin Smelt","Chinook Salmon", "Steelhead/Rainbow Trout")

#level_order = c("Green Sturgeon", "White Sturgeon","Delta Smelt", "Longfin Smelt", "Wakasagi","Chinook Salmon", "Steelhead/Rainbow Trout","Sacramento Splittail", "Sacramento Pikeminnow", "Sacramento Sucker", "Sacramento Blackfish", "Hitch", "Tule Perch", "Threespine Stickleback", "Prickly Sculpin", "Mississippi Silverside", "Striped Bass", "Largemouth Bass", "Giant Reed EAV", "Brazilian Waterweed", "Water Hyacinth FAV", "Microcystis", "Asian Clam", "Overbite Clam")
```

```{r}
data_exceedance <- data_exceedance0 %>%
  mutate(fWY = factor(WY),
         Status = ifelse(Species %in% native, "Native", "NonNative"),
         Status = as.factor(Status),
         Status2 = ifelse(Status == "Native", 0, 1),
         Lstatus = ifelse(Species %in% listed, "Listed", "Unlisted"),
         Lstatus = as.factor(Lstatus),
         Lstatus2 = ifelse(Lstatus == "Unlisted", 0, 1),
         daysTotal = ifelse(WY%%4 ==0, 366, 365))
```

# Adult days tolerance
## dataset
```{r}
data_adult_tol <- data_exceedance%>%
  mutate(log_daysTol = log(daysExceedanceTol + 1),
         sqrt_daysTol = sqrt(daysExceedanceTol)) %>%
  filter(LifeStage == "Adult" | LifeStage == "Undifferentiated", 
         !is.na(daysExceedanceTol)) %>%
  select(Status, Lstatus, fWY, Species, daysExceedanceTol, log_daysTol, sqrt_daysTol, daysTotal, Status2, Lstatus2) %>%
  unique()

str(data_adult_tol)
```

Write table of species included in model
```{r}
table <- data_exceedance %>% 
  filter(LifeStage == "Adult" | LifeStage == "Undifferentiated", 
         !is.na(daysExceedanceTol)) %>%
  select(Species, Lstatus, Tolerance_Upper) %>%
  distinct() %>%
  arrange(Species)

# write_csv(table, "SpeciesTables/data_write/model_status_species_thresholds.csv")
```



## normality
```{r}
hist(data_adult_tol$daysExceedanceTol)
hist(data_adult_tol$sqrt_daysTol)
hist(data_adult_tol$log_daysTol)
```

## plot categories
```{r}
ggplot(data_adult_tol, aes(x = Status, y = daysExceedanceTol)) + geom_boxplot() + theme(axis.text = element_text(size = 12))
ggplot(data_adult_tol, aes(x = Lstatus, y = daysExceedanceTol)) + geom_boxplot()

ggplot(data_adult_tol, aes(y = daysExceedanceTol)) + geom_histogram() + facet_wrap(~Lstatus, scales = "free")
ggplot(data_adult_tol, aes(y = sr_daysTol)) + geom_histogram() + facet_wrap(~Lstatus)
```

## t-test: significant
```{r}
t.test(data_adult_tol$Status2, data_adult_tol$sqrt_daysTol)
t.test(data_adult_tol$Lstatus2, data_adult_tol$daysExceedanceTol)
t.test(data_adult_tol$Lstatus2, data_adult_tol$sqrt_daysTol)
```

## kruskal wallis test
```{r}
kruskal.test(data_adult_tol$Status2, data_adult_tol$daysExceedanceTol)
kruskal.test(data_adult_tol$Lstatus2, data_adult_tol$daysExceedanceTol)
```


## lm
### by native status
```{r}
mAdult1 <- lm(log_daysTol ~ Status, data = data_adult_tol)
summary(mAdult1)
anova(mAdult1)
par(mfrow = c(2,2))
plot(mAdult1)
```

binomial
```{r}
mAdult1 <- glm(cbind(daysExceedanceTol, daysTotal) ~ Status, data = data_adult_tol, family = "binomial")
summary(mAdult1)
anova(mAdult1)
par(mfrow = c(2,2))
plot(mAdult1)

```



### by listed status
```{r}
mAdult2 <- lm(log_daysTol ~ Lstatus, data = data_adult_tol)
summary(mAdult2)
anova(mAdult2)
par(mfrow = c(2,2))
plot(mAdult2)
```

binomial
```{r}
mAdult2b <- glm(cbind(daysExceedanceTol, daysTotal) ~ Lstatus, data = data_adult_tol, family = binomial)
summary(mAdult2b)
anova(mAdult2b)
par(mfrow = c(2,2))
plot(mAdult2b)
```

## lmer

### by native status
```{r}
mAdultStat <- lmer(log_daysTol ~ Status + (1|Species), data = data_adult_tol)
summary(mAdultStat)
```

```{r}
mAdultStatusb <- glmer(cbind(daysExceedanceTol,daysTotal) ~ Status + (1|Species), data = data_adult_tol, family = "binomial")
summary(mAdultStatusb)

par(mfrow = c(2,3))
resid.tolb = resid(mAdultStatusb, type = "pearson")
lev2a = hatvalues(mAdultStatusb)
plot(mAdultStatusb)
qqnorm(resid(mAdultStatusb))
qqline(resid(mAdultStatusb))
plot(lev2a, y = resid.tolb)
hist(resid.tolb)
plot(resid.tolb)
plot(data_adult_tol$Status, resid.tolb)
acf(resid.tolb)

re_m_tolb <- as.data.frame(ranef(mAdultStatusb))

ggplot(re_m_tolb) + 
  geom_point(aes(grp, condval)) +
  geom_errorbar(aes(x = grp, ymin = condval-condsd, ymax = condval+condsd), width = 0.3) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 0.95),
      axis.title = element_blank())
```

### by listed status
```{r}
mAdultList <- lmer(sqrt_daysTol ~ Lstatus + (1|Species), data = data_adult_tol)
summary(mAdultList)
```

binomial
```{r}
mAdultListb <- glmer(cbind(daysExceedanceTol,daysTotal) ~ Lstatus + (1|Species), data = data_adult_tol, family = "binomial")
summary(mAdultListb)

par(mfrow = c(2,3))
resid.tolb = resid(mAdultListb, type = "pearson")
lev2a = hatvalues(mAdultListb)
plot(mAdultListb)
qqnorm(resid(mAdultListb))
qqline(resid(mAdultListb))
plot(lev2a, y = resid.tolb)
hist(resid.tolb)
plot(resid.tolb)
plot(data_adult_tol$Lstatus, resid.tolb)
acf(resid.tolb)

re_m_tolb <- as.data.frame(ranef(mAdultListb))

ggplot(re_m_tolb) + 
  geom_point(aes(grp, condval)) +
  geom_errorbar(aes(x = grp, ymin = condval-condsd, ymax = condval+condsd), width = 0.3) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 0.95),
      axis.title = element_blank())

```

```{r}
library(effects)
library(sjPlot)

c3 <- with(effect("Lstatus",mAdultListb),cbind(lower,upper))

(mem_eff=Effect(focal.predictors = "Lstatus", mod=mAdultListb))
eff.lstatus <- allEffects(mod = mAdultListb, fixed.predictors = "Lstatus")
eff.lstatus
tab_model(mAdultListb)
confint(mAdultListb)
```

#### post hoc
```{r}
mtol.emm <- emmeans(mAdultListb , "Lstatus", data=data_adult_tol)
pairs(mtol.emm, adjust = "tukey")
plot(mtol.emm, comparisons = TRUE, type = "response") + theme_bw()
confint(mtol.emm, type = "response")
```


# Adult days optimum
## dataset
```{r}
data_adult_opt <- data_exceedance %>%
  mutate(log_daysOpt = log(daysExceedanceOpt + 1),
         sqrt_daysOpt = sqrt(daysExceedanceOpt)) %>%
  filter(LifeStage == "Adult" | LifeStage == "Undifferentiated",  
         !is.na(daysExceedanceOpt)) %>%
  select(Status, Lstatus, Status2, Lstatus2, fWY, Species, daysExceedanceOpt, log_daysOpt, sqrt_daysOpt, daysTotal) %>%
  unique() 

str(data_adult_opt)
```

Write table of species included in model
```{r}
table_opt <- data_exceedance %>% 
  filter(LifeStage == "Adult" | LifeStage == "Undifferentiated", 
         !is.na(daysExceedanceOpt)) %>%
  select(Species, Lstatus, Suboptimum_Upper) %>%
  distinct() %>%
  arrange(Species)

# write_csv(table_opt, "SpeciesTables/data_write/model_status_species_thresholds_opt.csv")
```



## normality
```{r}
par(mfrow = c(1,1))
hist(data_adult_opt$daysExceedanceOpt)
hist(data_adult_opt$sqrt_daysOpt)
hist(data_adult_opt$log_daysOpt)
```

## plot categories
```{r}
ggplot(data_adult_opt, aes(x = Status, y = daysExceedanceOpt)) + geom_boxplot() + theme(axis.text = element_text(size = 12))
ggplot(data_adult_opt, aes(x = Lstatus, y = daysExceedanceOpt)) + geom_boxplot()+ theme(axis.text = element_text(size = 12))

ggplot(data_adult_opt, aes(y = daysExceedanceOpt)) + geom_histogram() + facet_wrap(~Lstatus)
```

## t-test: significant
```{r}
t.test(data_adult_opt$Status2, data_adult_opt$daysExceedanceOpt)
t.test(data_adult_opt$Lstatus2, data_adult_opt$daysExceedanceOpt)
```

## kruskal wallis: significant
```{r}
kruskal.test(data_adult_opt$Status2, data_adult_opt$daysExceedanceOpt)
kruskal.test(data_adult_opt$Lstatus2, data_adult_opt$daysExceedanceOpt)
```


## lm
### by native status
```{r}
mAdult1_o <- lm(daysExceedanceOpt ~ Status, data = data_adult_opt)
summary(mAdult1_o)
anova(mAdult1_o)
par(mfrow = c(2,2))
plot(mAdult1_o)
```

```{r}
mAdult1b_o <- glm(cbind(daysExceedanceOpt, daysTotal) ~ Status, data = data_adult_opt, family = binomial)
summary(mAdult1b_o)
anova(mAdult1b_o)
par(mfrow = c(2,2))
plot(mAdult1b_o)
```


### by listed status
```{r}
mAdult2_o <- lm(daysExceedanceOpt ~ Lstatus, data = data_adult_opt)
summary(mAdult2_o)
anova(mAdult2_o)
par(mfrow = c(2,2))
plot(mAdult2_o)
```

```{r}
mAdult2b_o <- glm(cbind(daysExceedanceOpt, daysTotal) ~ Lstatus, data = data_adult_opt, family = binomial)
summary(mAdult2b_o)
anova(mAdult2b_o)
par(mfrow = c(2,2))
plot(mAdult2b_o)
```

## lmer

### by native status
```{r}
mAdultStat_o <- lmer(daysExceedanceOpt ~ Status + (1|Species), data = data_adult_opt)
summary(mAdultStat_o)


par(mfrow = c(2,3))
resid.opt = resid(mAdultStat_o, type = "pearson")
lev2a = hatvalues(mAdultStat_o)
plot(mAdultStat_o)
qqnorm(resid(mAdultStat_o))
qqline(resid(mAdultStat_o))
plot(lev2a, y = resid.opt)
hist(resid.opt)
plot(resid.opt)
plot(data_adult_opt$Status, resid.opt)
acf(resid.opt)
```

binomial
```{r}
mAdultStat_ob <- glmer(cbind(daysExceedanceOpt,daysTotal) ~ Status +  (1|Species), data = data_adult_opt, family = "binomial")
summary(mAdultStat_ob)

par(mfrow = c(2,3))
resid.optb = resid(mAdultStat_ob, type = "pearson")
lev2a = hatvalues(mAdultStat_ob)
plot(mAdultStat_ob)
qqnorm(resid(mAdultStat_ob))
qqline(resid(mAdultStat_ob))
plot(lev2a, y = resid.optb)
hist(resid.optb)
plot(resid.optb)
plot(data_adult_opt$Status, resid.optb)
acf(resid.optb)

re_m_ob <- as.data.frame(ranef(mAdultStat_ob))

ggplot(re_m_ob) + 
  geom_point(aes(grp, condval)) +
  geom_errorbar(aes(x = grp, ymin = condval-condsd, ymax = condval+condsd), width = 0.3) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 0.95),
      axis.title = element_blank())
```


### by listed status
```{r}
mAdultList_o <- lmer(daysExceedanceOpt ~ Lstatus + (1|Species), data = data_adult_opt)
summary(mAdultList_o)

par(mfrow = c(2,3))
resid.opt = resid(mAdultList_o, type = "pearson")
lev2a = hatvalues(mAdultList_o)
plot(mAdultList_o)
qqnorm(resid(mAdultList_o))
qqline(resid(mAdultList_o))
plot(lev2a, y = resid.opt)
hist(resid.opt)
plot(resid.opt)
plot(data_adult_opt$Lstatus, resid.opt)
acf(resid.opt)
```

binomial
```{r}
mAdultList_ob <- glmer(cbind(daysExceedanceOpt,daysTotal) ~ Lstatus +  (1|Species), data = data_adult_opt, family = "binomial")
summary(mAdultList_ob)

summary(mAdult2b_o)




par(mfrow = c(2,3))
resid.optb = resid(mAdultList_ob, type = "pearson")
lev2a = hatvalues(mAdultList_ob)
plot(mAdultList_ob)
qqnorm(resid(mAdultList_ob))
qqline(resid(mAdultList_ob))
plot(lev2a, y = resid.optb)
hist(resid.optb)
plot(resid.optb)
plot(data_adult_opt$Lstatus, resid.optb)
acf(resid.optb)

re_m_ob <- as.data.frame(ranef(mAdultList_ob))

ggplot(re_m_ob) + 
  geom_point(aes(grp, condval)) +
  geom_errorbar(aes(x = grp, ymin = condval-condsd, ymax = condval+condsd), width = 0.3) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 0.95),
      axis.title = element_blank())
```

#### post hoc
```{r}
mopt.emm <- emmeans(mAdultList_ob , "Lstatus", data=data_adult_opt)
pairs(mopt.emm, adjust = "tukey")
plot(mopt.emm, comparisons = TRUE, type = "response") + theme_bw()
confint(mopt.emm, type = "response")
```


# Tolerance Values
```{r}
tol <- data_exceedance %>%
  filter(LifeStage == "Adult" | LifeStage == "Undifferentiated", 
         !is.na(Tolerance_Upper)) %>%
  select(Status,Status2,  Lstatus, Lstatus2, Species, Tolerance_Upper) %>%
  unique()
```

Not normally distributed...
```{r}
par(mfrow = c(1,1))
hist(tol$Tolerance_Upper)
hist(sqrt(tol$Tolerance_Upper))
```

## t-test
```{r}
t.test(tol$Status2, tol$Tolerance_Upper)
t.test(tol$Lstatus2, tol$Tolerance_Upper)
```

## kruskal wallis: ns
```{r}
kruskal.test(tol$Status2, tol$Tolerance_Upper)
kruskal.test(tol$Lstatus2, tol$Tolerance_Upper)
```

## lm

### Native status: non-significant
```{r}
mTol <- lm(Tolerance_Upper ~ Status, data = tol)
summary(mTol)
anova(mTol)
par(mfrow = c(2,2))
plot(mTol)
```

### Listed status: significant
```{r}
mTol2 <- lm(Tolerance_Upper ~ Lstatus, data = tol)
summary(mTol2)
anova(mTol2)
par(mfrow = c(2,2))
plot(mTol2)
par(mfrow = c(1,1))
hist(resid(mTol2, "pearson"))
```


# Optimum Values - only 17 values
```{r}
opt <- data_exceedance %>%
  filter(LifeStage == "Adult" | LifeStage == "Undifferentiated", 
         !is.na(Suboptimum_Upper)) %>%
  select(Status,Status2,  Lstatus, Lstatus2, Species, Suboptimum_Upper) %>%
  unique()
```

Not normally distributed...
```{r}
hist(opt$Suboptimum_Upper)
hist(sqrt(opt$Suboptimum_Upper))
```

## t-test: Significant
```{r}
t.test(opt$Status2, opt$Suboptimum_Upper)
t.test(opt$Lstatus2, opt$Suboptimum_Upper)
```

## kruskal wallis: ns
```{r}
kruskal.test(opt$Status2, opt$Suboptimum_Upper)
kruskal.test(opt$Lstatus2, opt$Suboptimum_Upper)
```

## lm
### Native status: significant
```{r}
mOpt <- lm(Suboptimum_Upper ~ Status, data = opt)
summary(mOpt)
anova(mOpt)
par(mfrow = c(2,2))
plot(mOpt)
hist(resid(mOpt, "pearson"))
```

### Listed status: significant
```{r}
mOpt2 <- lm(Suboptimum_Upper ~ Lstatus, data = opt)
summary(mOpt2)
anova(mOpt2)
par(mfrow = c(2,2))
plot(mOpt2)
par(mfrow = c(1,1))
hist(resid(mOpt2, "pearson"))
```