lme_alpha.Rmd

---
title: "Modeling Changes in Alpha Diversity over Time"
output:
  workflowr::wflow_html:
    toc: true
editor_options:
  chunk_output_type: console
---

This page contains the investigation of the changes in alpha diversity metrics (Observed OTUs, Shannon, and Inverse Simpson) over time.

```{r data, include=FALSE, echo=FALSE}
# load packages
source("code/load_packages.R")

# get data
source("code/get_cleaned_data.R")

theme_set(theme_bw())
pal = "Set1"
scale_colour_discrete <-  function(palname=pal, ...){
  scale_colour_brewer(palette=palname, ...)
}
scale_fill_discrete <-  function(palname=pal, ...){
  scale_fill_brewer(palette=palname, ...)
}

knitr::opts_chunk$set(out.width = "225%")

```


Due to only having 11 people, we will only estimate 1 random effect per model (random intercepts).
All other effects are fixed, but we will only estimate up to 2 effects.
Meaning we can only include two covariates in each model at a time.

Here, we investigated how the metrics of alpha diversity changed over time.
The general model tested for which will have reported results in the manuscript is a growth model of each alpha diversity metric.
The final growth model had the following form:

\begin{align*}
Y_{ij}&=\beta_{0j} + \beta_{1j}(Week)_{ij} + r_{ij}\\
\beta_{0j}&= \gamma_{00} + \gamma_{01}(Intervention)_{j} + u_{0j}\\
\beta_{1j}&= \gamma_{10} + \gamma_{11}(Intervention)_{j}\\
\end{align*}
where, the major aim was to assess the significance of two terms:

1. the main effect of the intervention on the intercept - indicating average changes in alpha diversity between intervention groups.
2. the interaction between time and intervention (intB:time), which assessing the magnitude of the change in alpha diversity as the intervention progresses. 

```{r gen-code}

ICC <- function(x){
  icc <- VarCorr(x)[[1]]/(VarCorr(x)[[1]] + sigma(x)**2)
  icc <- lapply(icc, function(x) { attributes(x) <- NULL; x })
  icc <- icc[[1]]
  return(icc)
}

mydata <- microbiome_data$meta.dat %>%
  mutate(intB = ifelse(Intervention=="B", 1,0),
         time = as.numeric(Week) - 1,
         female = ifelse(Gender == "F", 1, 0),
         hispanic = ifelse(Ethnicity %in% c("White", "Asian", "Native America"), 1, 0))

# Objects to say results
alpha.fit <- list()

```


## Alpha: Observed OTUs

### Unconditional Model

```{r}

# lmer - for alpha metrics
# unconditional model
fit <- lmer(Observed ~ 1 + (1 | SubjectID),
            data = mydata)
summary(fit)
plot(fit)
ICC(fit)

```

### Fixed Effect of Time

```{r}


fit <- lmer(Observed ~ 1 + time + (1 | SubjectID),
            data = mydata)
summary(fit)
plot(fit)
# plot
dat <- cbind(mydata, fit=predict(fit))
ggplot(dat, aes(time, Observed, group=SubjectID))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```

### Fixed and Random effect of Time

```{r}


fit <- lmer(Observed ~ 1 + time + (1 + time || SubjectID),
            data = mydata)
summary(fit)
ranova(fit) # can take out random effect time
plot(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Observed, group=SubjectID))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)#+
  #geom_abline(intercept = fixef(fit)[1], slope=fixef(fit)[2],
              #linetype="dashed", size=1.5)

```


### Adding demographics

Time is only a fixed effect.

```{r}


fit <- lmer(Observed ~ 1 + time + female + hispanic + (1 | SubjectID),
            data = mydata)
summary(fit)
ranova(fit)
plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Observed, group=SubjectID, color=Gender))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```


### Intervention effect

```{r}


fit <- lmer(Shannon ~ 1 + intB + time + female +  (1 | SubjectID),
            data = mydata)
summary(fit)
ranova(fit)
plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Shannon, group=SubjectID, color=Intervention))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```


#### Overall Effect

```{r}


fit <- lmer(Observed ~ 1 + intB + female + hispanic + (1 | SubjectID),
            data = mydata)
summary(fit)

plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Observed, group=SubjectID, color=Intervention))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```

#### With Time

```{r}

fit <- lmer(Observed ~ 1 + intB*time + (1 | SubjectID),
            data = mydata)
alpha.fit[["Observed"]]<- summary(fit)[["coefficients"]]
summary(fit)
ranova(fit)
plot(fit)
ICC(fit)


dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Observed, group=SubjectID, color=Intervention))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)
```


## Alpha: Shannon Diversity

### Unconditional Model

```{r}

# lmer - for alpha metrics
# unconditional model
fit <- lmer(Shannon ~ 1 + (1 | SubjectID),
            data = mydata)
summary(fit)
plot(fit)
ICC(fit)

```

### Fixed Effect of Time

```{r}


fit <- lmer(Shannon ~ 1 + time + (1 | SubjectID),
            data = mydata)
summary(fit)
plot(fit)
# plot
dat <- cbind(mydata, fit=predict(fit))
ggplot(dat, aes(time, Shannon, group=SubjectID))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```

### Fixed and Random effect of Time

```{r}


fit <- lmer(Shannon ~ 1 + time + (1 + time || SubjectID),
            data = mydata)
summary(fit)
ranova(fit) # can take out random effect time
plot(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Shannon, group=SubjectID))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)#+
  #geom_abline(intercept = fixef(fit)[1], slope=fixef(fit)[2],
              #linetype="dashed", size=1.5)

```


### Adding demographics

Time is only a fixed effect.

```{r}


fit <- lmer(Shannon ~ 1 + time + female + hispanic + (1 | SubjectID),
            data = mydata)
summary(fit)
ranova(fit)
plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Shannon, group=SubjectID, color=Gender))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```


### Intervention effect


```{r}


fit <- lmer(Shannon ~ 1 + intB + time + female +  (1 | SubjectID),
            data = mydata)
summary(fit)
ranova(fit)
plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Shannon, group=SubjectID, color=Intervention))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```

#### Overall Effect

```{r}


fit <- lmer(Shannon ~ 1 + intB + female + hispanic + (1 | SubjectID),
            data = mydata)
summary(fit)

plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Shannon, group=SubjectID, color=Intervention))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```

#### With Time

```{r}

fit <- lmer(Shannon ~ 1 + intB*time + (1 | SubjectID),
            data = mydata)
alpha.fit[["Shannon"]]<- summary(fit)[["coefficients"]]
summary(fit)
ranova(fit)
plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, Shannon, group=SubjectID, color=Intervention))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)
```



## Alpha: Inverse Simpson

### Unconditional Model

```{r}

# lmer - for alpha metrics
# unconditional model
fit <- lmer(InvSimpson ~ 1 + (1 | SubjectID),
            data = mydata)
summary(fit)
plot(fit)
ICC(fit)

```

### Fixed Effect of Time

```{r}


fit <- lmer(InvSimpson ~ 1 + time + (1 | SubjectID),
            data = mydata)
summary(fit)
plot(fit)
# plot
dat <- cbind(mydata, fit=predict(fit))
ggplot(dat, aes(time, InvSimpson, group=SubjectID))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```

### Fixed and Random effect of Time

```{r}


fit <- lmer(InvSimpson ~ 1 + time + (1 + time || SubjectID),
            data = mydata)
summary(fit)
ranova(fit) # can take out random effect time
plot(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, InvSimpson, group=SubjectID))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)#+
  #geom_abline(intercept = fixef(fit)[1], slope=fixef(fit)[2],
              #linetype="dashed", size=1.5)

```


### Adding demographics

Time is only a fixed effect.

```{r}


fit <- lmer(InvSimpson ~ 1 + time + female + hispanic + (1 | SubjectID),
            data = mydata)
summary(fit)
ranova(fit)
plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, InvSimpson, group=SubjectID, color=Gender))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```


### Intervention effect


```{r}


fit <- lmer(InvSimpson ~ 1 + intB + time + female +  (1 | SubjectID),
            data = mydata)
summary(fit)
ranova(fit)
plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, InvSimpson, group=SubjectID, color=Intervention))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```


#### Overall Effect

```{r}


fit <- lmer(InvSimpson ~ 1 + intB + female + hispanic + (1 | SubjectID),
            data = mydata)
summary(fit)

plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, InvSimpson, group=SubjectID, color=Intervention))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)

```

#### With Time

```{r}

fit <- lmer(InvSimpson ~ 1 + intB*time + (1 | SubjectID),
            data = mydata)
alpha.fit[["InvSimpson"]]<- summary(fit)[["coefficients"]]
summary(fit)
ranova(fit)
plot(fit)
ICC(fit)

dat <- cbind(mydata, fit=predict(fit))

ggplot(dat, aes(time, InvSimpson, group=SubjectID, color=Intervention))+
  geom_line(aes(y=fit))+
  geom_point(alpha=0.5)
```



## Summary of Results

To recap, the final model tested for which will have reported results in the manuscript is a growth model of each alpha diversity metric.
The final growth model had the following form:

\begin{align*}
Y_{ij}&=\beta_{0j} + \beta_{1j}(Week)_{ij} + r_{ij}\\
\beta_{0j}&= \gamma_{00} + \gamma_{01}(Intervention)_{j} + u_{0j}\\
\beta_{1j}&= \gamma_{10} + \gamma_{11}(Intervention)_{j}\\
\end{align*}
where, the major aim was to assess the significance of two terms:

1. the main effect of the intervention on the intercept - indicating average changes in alpha diversity between intervention groups.
2. the interaction between time and intervention (intB:time), which assessing the magnitude of the change in alpha diversity as the intervention progresses. 

```{r summary, echo=F}

unlist_results <- function(fit){
  out <- fit[[1]]
  results.out <- data.frame(matrix(nrow=nrow(out), ncol=ncol(out)+2))
  colnames(results.out) <- c("Metric", "Term", colnames(out))
  results.out[,1] <- names(fit)[1]
  results.out[,2] <- rownames(out)
  results.out[,3:7] <- out[,1:5]
  
  i <-  2
  for(i in 2:length(fit)){
    out <- fit[[i]]
    results.i <- data.frame(matrix(nrow=nrow(out), ncol=ncol(out)+2))
    colnames(results.i) <- c("Metric", "Term", colnames(out))
    results.i[,1] <- names(fit)[i]
    results.i[,2] <- rownames(out)
    results.i[,3:7] <- out[,1:5]
    results.out <- full_join(results.out, results.i)
  }
  # exclude intercept
  results.out <- filter(results.out, Term != "(Intercept)")
  
  results.out$p.adjust <- p.adjust(results.out$`Pr(>|t|)`, "fdr")
  
  return(results.out)
}
```

```{r}

alpha.results <- unlist_results(alpha.fit)

kable(alpha.results, format="html", digits=3)%>%
  kable_styling(full_width = T)

```

Therefore, the intervention did not have a statistically significant effect on alpha diversity metrics.