---
title: "Sham/Ovex Time Series Analyis"
author: "Christine Olmstead"
date: "4/12/2020"
output: html_document
---

```{r setup}
    library(dplyr)
    library(ggplot2)
    library(ggfortify)
    library(lme4)
    library(nlme)
    library(car)
    library(MuMIn)
knitr::opts_chunk$set(echo = TRUE)
```

## Faiths_PD

```{r}
#---------------------------faiths_Pd-----------------------------------
ovex.pd <- read.csv("faith_pd_sham_ovex.csv")
  
ovex.filtered.pd <- ovex.pd %>% filter(time != 0)

#Filter by time 
ovex.T0.pd <- ovex.pd %>% filter(time == 0)
ovex.T1.pd <- ovex.pd %>% filter(time == 1)
ovex.T2.pd <- ovex.pd %>% filter(time == 2)
ovex.T3.pd <- ovex.pd %>% filter(time == 3)
ovex.T4.pd <- ovex.pd %>% filter(time == 4)
ovex.T5.pd <- ovex.pd %>% filter(time == 5)

#Filter by treatment 
ovex.o.pd <- ovex.pd %>% filter(Treatment == "Ovex")
ovex.s.pd <- ovex.pd %>% filter(Treatment == "Sham")

#Check for normality
ks.test(ovex.pd$faith_pd, "pnorm", mean=mean(ovex.pd$faith_pd), sd=sd(ovex.pd$faith_pd))

shapiro.test(ovex.pd$faith_pd)

#--------------Plot Sham Vs Ovex by Time-----------------------
ggplot(ovex.pd, aes(x = Treatment, y=faith_pd, fill = Treatment)) +
  geom_boxplot(aes(middle = mean(faith_pd))) +
  scale_fill_manual(values=c("grey", "white")) +
  facet_wrap(~time) +
 theme_light() + 
  theme(axis.line.x = element_line(color="black", size = 0.7),axis.text.x = element_text(size=12, family="Arial"),axis.text.y = element_text(size=12,  family="Arial"),text = element_text(size=14, family="Arial"),panel.grid.major = element_blank(),panel.grid.minor = element_blank(), plot.margin=unit(c(.5,.5,.5,.5),"cm"),  plot.title = element_text(hjust = 0.2, vjust = 2)) + labs(title = "Faith's Phylogenetic Diversity (week 0-5)", y = "Faith's PD", x = "")


#-------------------------Plot Ovex over time -----------------------------------
ggplot(ovex.o.pd, aes(x = as.factor(time), y=faith_pd)) +
  geom_boxplot() +
  ggtitle("Faiths PD vs. Time Ovex") +
  xlab("Time") + 
  ylab("Faiths PD") +
  theme_bw()

#---------------------------Plot Sham over time------------------------------------
ggplot(ovex.s.pd, aes(x = as.factor(time), y=faith_pd)) +
  geom_boxplot() +
  ggtitle("Faiths PD vs. Time Sham") +
  xlab("Time") + 
  ylab("Faiths PD") +
  theme_bw()

#---------------------Modeling---------------------------------------------

#Compare ovex and sham for ecah timepoint 
model.t0 <- lm(faith_pd ~ Treatment, data= ovex.T0.pd, na.action = na.omit)
anova(model.t0)
summary(model.t0)

model.t1 <- lm(faith_pd ~ Treatment, data= ovex.T1.pd, na.action = na.omit)
anova(model.t1)
summary(model.t1)

model.t2 <- lm(faith_pd ~ Treatment, data= ovex.T2.pd, na.action = na.omit)
anova(model.t2)
summary(model.t2)

model.t3 <- lm(faith_pd ~ Treatment, data= ovex.T3.pd, na.action = na.omit)
anova(model.t3)
summary(model.t3)

model.t4 <- lm(faith_pd ~ Treatment, data= ovex.T4.pd, na.action = na.omit)
anova(model.t4)
summary(model.t4)

model.t5 <- lm(faith_pd ~ Treatment, data= ovex.T5.pd, na.action = na.omit)
anova(model.t5)
summary(model.t5)

#Ovex modeling 
ggplot(ovex.o.pd, aes(faith_pd)) + 
  geom_histogram(binwidth = 1) +
  xlab( "Faith's PD Ovex") +
  ylab("Count") +
  theme_bw()

model1 <- lme(faith_pd ~ time, random= ~1| iden, data= ovex.o.pd, method = "ML", na.action = na.omit)

plot(model1)
qqnorm(model1)

anova(model1)
summary(model1)

#Sham modeling
ggplot(ovex.s.pd, aes(faith_pd)) + 
  geom_histogram(binwidth = 1) +
  xlab( "Faith's PD Sham") +
  ylab("Count") +
  theme_bw()

model2 <- lme(faith_pd ~ time, random= ~1| iden, data= ovex.s.pd, method = "ML", na.action = na.omit)

plot(model2)
qqnorm(model2)

anova(model2)


#model and include treatment and time 
model3 <- nlme::lme(faith_pd ~ Treatment + time, random= ~1| iden, data= ovex.pd, method = "ML", na.action = na.omit)
anova(model3, test = "chi")
summary(aov(model3))
summary(model3)
r.squaredGLMM(model3)



```

## Observed SVs


```{r}

ovex.sv <- read.csv("observed_svs_sham_ovex.csv")
  
ovex.filtered.sv <- ovex.sv %>% filter(time != 0)

#filter by time point 
ovex.T0.sv <- ovex.sv %>% filter(time == 0)
ovex.T1.sv <- ovex.sv %>% filter(time == 1)
ovex.T2.sv <- ovex.sv %>% filter(time == 2)
ovex.T3.sv <- ovex.sv %>% filter(time == 3)
ovex.T4.sv <- ovex.sv %>% filter(time == 4)
ovex.T5.sv <- ovex.sv %>% filter(time == 5)

#Filter by treatment 
ovex.o.sv <- ovex.sv %>% filter(Treatment == "Ovex")
ovex.s.sv <- ovex.sv %>% filter(Treatment == "Sham")

shapiro.test(ovex.sv$observed_otus)

#--------------Plot Sham Vs Ovex by Time-----------------------
ggplot(ovex.sv, aes(x = Treatment, y=observed_otus, fill = Treatment)) +
  geom_boxplot(aes(middle = mean(observed_otus))) +
  scale_fill_manual(values=c("grey", "white")) +
  ggtitle("Observed SVs vs. Treatment") +
  xlab("Treatment") + 
  ylab("Observed SVs") +
  facet_wrap(~time) +
  theme_bw()

ggplot(ovex.sv, aes(x = Treatment, y=observed_otus, fill = Treatment)) +
  geom_boxplot(aes(middle = mean(observed_otus))) +
  scale_fill_manual(values=c("grey", "white")) +
  facet_wrap(~time) +
 theme_light() + 
  theme(axis.line.x = element_line(color="black", size = 0.7),axis.text.x = element_text(size=12, family="Arial"),axis.text.y = element_text(size=12,  family="Arial"),text = element_text(size=14, family="Arial"),panel.grid.major = element_blank(),panel.grid.minor = element_blank(), plot.margin=unit(c(.5,.5,.5,.5),"cm"),  plot.title = element_text(hjust = 0.2, vjust = 2)) + labs(title = "Observed Sequence Variants (week 0-5)", y = "Observed SVs", x = "")


#-------------------------Plot Ovex over time -----------------------------------
ggplot(ovex.o.sv, aes(x = as.factor(time), y=observed_otus)) +
  geom_boxplot() +
  ggtitle("Observed SVs vs. Time Ovex") +
  xlab("Time") + 
  ylab("Observed SVs") +
  theme_bw()

#-------------------------Plot Sham over time -----------------------------------
ggplot(ovex.s.sv, aes(x = as.factor(time), y=observed_otus)) +
  geom_boxplot() +
  ggtitle("Observed SVs vs. Time Sham") +
  xlab("Time") + 
  ylab("Observed SVs") +
  theme_bw()

#---------------------Modeling---------------------------------------------

#Compare ovex and sham for ecah timepoint 
model.sv.t0 <- lm(observed_otus ~ Treatment, data= ovex.T0.sv, na.action = na.omit)
anova(model.sv.t0)
summary(model.sv.t0)

model.sv.t1 <- lm(observed_otus ~ Treatment, data= ovex.T1.sv, na.action = na.omit)
anova(model.sv.t1)
summary(model.sv.t1)

model.sv.t2 <- lm(observed_otus ~ Treatment, data= ovex.T2.sv, na.action = na.omit)
anova(model.sv.t2)
summary(model.sv.t2)

model.sv.t3 <- lm(observed_otus ~ Treatment, data= ovex.T0.sv, na.action = na.omit)
anova(model.sv.t3)
summary(model.sv.t3)

model.sv.t4 <- lm(observed_otus ~ Treatment, data= ovex.T4.sv, na.action = na.omit)
anova(model.sv.t4)
summary(model.sv.t4)

model.sv.t5 <- lm(observed_otus ~ Treatment, data= ovex.T5.sv, na.action = na.omit)
anova(model.sv.t5)
summary(model.sv.t5)



#Ovex modeling 
ggplot(ovex.o.sv, aes(observed_otus)) + 
  geom_histogram(binwidth = 10) +
  xlab( "Observed SVs Ovex") +
  ylab("Count") +
  theme_bw()

model1.sv <- lme(observed_otus ~ as.factor(time), random= ~1| iden, data= ovex.o.sv, method = "ML", na.action = na.omit)

plot(model1.sv)
qqnorm(model1.sv)

anova(model1.sv)

summary(model1.sv)



#Sham modeling 
ggplot(ovex.s.sv, aes(observed_otus)) + 
  geom_histogram(binwidth = 10) +
  xlab( "Observed SV Sham") +
  ylab("Count") +
  theme_bw()

model2.sv <- lme(observed_otus ~ as.factor(time), random= ~1| iden, data= ovex.s.sv, method = "ML", na.action = na.omit)

plot(model2.sv)
qqnorm(model2.sv)

anova(model2.sv)

#model and include treatment and time 
model3.sv <- lme(observed_otus ~ Treatment + time, random= ~1| iden, data= ovex.sv, method = "ML", na.action = na.omit)

anova(model3.sv)
summary(model3.sv)
r.squaredGLMM(model3.sv)

```


