---
title: "Model presence of MOTUs in LOWA diet with GLMM"
author: "Nathan Brouwer"
date: "February 23, 2018"
output: html_document
editor_options: 
  chunk_output_type: console
---

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

## Introduction
This script uses a GLMM to model the relationship between the presence of MOTUs in LOWA nestling diets and EPT availability.

## Preliminaries

### Load libraries

```{r}
library(lme4)
library(arm) #for se.coef() function to get se for blups

library(merTools) #for predictInterval() for plotting

library(stringr)#for cleaning MOTU names

library(ggplot2)
library(cowplot)


```

### Load data

```{r}
dat_mult_spp_glmm <- read.csv(file = "data_mult_spp_glmm.csv")
```

### Set up predictor

Extract genera names.  

```{r}
#split up genera column
family. <- str_split(dat_mult_spp_glmm$family,"\\.",
                     simplify = T)

#create columns for each taxonomic level
dat_mult_spp_glmm$Class = family.[,2]
dat_mult_spp_glmm$Order  = family.[,3]
dat_mult_spp_glmm$Family = family.[,4]

#create factor varible for EPT
EPT. <-c("Ephemeroptera",
         "Plecoptera",
         "Trichoptera")

dat_mult_spp_glmm$EPT <- ifelse(dat_mult_spp_glmm$Order %in% EPT.,
                        "EPT",
                        "non-EPT")
```


## Run GLMM

* Note (1+ PER_EPT_cent|genera) for seperate slopes for each
genera (MOTU)
* (1|SAMPLE) is random effect for each sample since many genera occur w/in a single sample

### Pure random effects models
```{r}
# Null model - no random intercepts
m.null <- glmer(pres.abs ~ 1 + 
             (1|SAMPLE) +
             (1|STREAM) + 
             (1|NEST) + 
             (1|NEST:VISIT) +
             (1|family),
           data = dat_mult_spp_glmm,
           control = glmerControl(optimizer = "Nelder_Mead"),
           family = binomial)

# Random intercept for %EPT
m.per.ept <- glmer(pres.abs ~ 1 + 
             (1|SAMPLE) +
             (1|STREAM) + 
             (1|NEST) + 
             (1|NEST:VISIT) +
             (1 + PROP_EPT_cent|family),
           data = dat_mult_spp_glmm,
           control = glmerControl(optimizer = "Nelder_Mead"),
           family = binomial)
```



## Extract genera-specific slopes

* aka "blups" - best linear unbiased estiamte
* get SEs using arm::se.coef()
```{r}
blups.per.ept <- data.frame(slope = ranef(m.per.ept)$family$PROP_EPT_cent,
           SE    = se.coef(m.per.ept)[["family"]][, "PROP_EPT_cent"])

```


### Clean output

Extract genera names.  SHould be idential in the models but will do sep just in case
```{r}
family.ept <- str_split(rownames(blups.per.ept),"\\.",simplify = T)
rownames(blups.per.ept) <- NULL
```



### Build up df of BLUPS
```{r}
blups.per.ept <- data.frame( 
           Class = family.ept[,2],
           Order  = family.ept[,3],
           Family = family.ept[,4],
           blups.per.ept)
```



### Code EPT vs non-ept
```{r}
blups.per.ept$EPT <- ifelse(blups.per.ept$Order %in% c("Ephemeroptera",
                                               "Plecoptera",
                                               "Trichoptera"),
                        "EPT","non-EPT")
```




### +/- 1 SE

```{r}
blups.per.ept$SE.lo <- with(blups.per.ept, slope - 1*SE)
blups.per.ept$SE.hi <- with(blups.per.ept, slope + 1*SE)
```


### Approximate CIs using 1.96*SE
```{r}
blups.per.ept$CI.lo <- with(blups.per.ept, slope - 1.96*SE)
blups.per.ept$CI.hi <- with(blups.per.ept, slope + 1.96*SE)
```


### Approximate "significance"
Greater than 1 CI
```{r}
blups.per.ept$sig <- ifelse(blups.per.ept$slope > 0 & blups.per.ept$CI.lo > 0 |
       blups.per.ept$slope < 0 & blups.per.ept$CI.hi < 0, "*","" )
```


### W/in 1 SE
```{r}
blups.per.ept$mar.sig <- ifelse(blups.per.ept$slope > 0 & blups.per.ept$SE.lo > 0 |
       blups.per.ept$slope < 0 & blups.per.ept$SE.hi < 0, "*","" )
```


## Plot

### Set up plot

 Make sep copy in case I screw it up
```{r}
blups.per.ept2 <- blups.per.ept
```

##### Order by taxa

order 1st by Order, then Family
```{r}
#ept model
lvs.ept <- blups.per.ept2$Family[order(blups.per.ept2$Order,
                                       blups.per.ept2$Family)]
blups.per.ept2$Family <- factor(blups.per.ept2$Family,
                            levels = lvs.ept)

```



#### Select what to plot
 . >95% CI or > 1 SE from 0
```{r}
i.ept.sig <- which(blups.per.ept2$sig == "*")

i.ept.marg.sig <- which(blups.per.ept2$mar.sig == "*")



```

### Plot

#### Plot Percent EPT
```{r}
gg.ept <- ggplot(data = blups.per.ept2[i.ept.marg.sig,],
       aes(y = slope,
           x = Family,
           color = EPT,
           shape = Order)) +
  coord_flip() +
  geom_errorbar(aes(ymin = CI.lo,
                    ymax = CI.hi),
                width = 0) +
  geom_errorbar(aes(ymin = SE.lo,
                    ymax = SE.hi),
                width = 0,
                size = 2) +
  geom_point(size = 4) +
  scale_shape_manual(values=c(0,1,2,5,6,7,8,9,10,15,16,17)) +
  geom_hline(yintercept = 0) +
  theme(axis.text.y = element_text(size = 8)) +
  ylab("Effect Size (log odds ratio)") +
  xlab("MOTU Family") +
  geom_point(size = 2,shape = 20)

```


### Save plot
```{r}
ggsave(gg.ept,filename = "MOTUs_vs_EPT.png",
       width = 5,height = 8)
```