---
title: "R codes for: Diversity loss from multiple interacting disturbances is regime-dependent."
date: "Updated: `r format(Sys.time(), '%B %d, %Y')`"
author: "L. Lear, H. Inamine, K. Shea, and A. Buckling"
header-includes:
  - \newcommand{\freq}{\texttt{pulse frequency}}
  - \newcommand{\Freq}{\texttt{Pulse frequency}}
  - \newcommand{\inv}{\texttt{invasion}}
  - \newcommand{\Inv}{\texttt{Invasion}}
  - \newcommand{\intx}{\freq{}$\times$\inv{}}
  - \newcommand{\Intx}{\Freq{}$\times$\inv{}}
  - \newcommand{\dsq}{\texttt{I(disturbance\string^2)}}
output:
  pdf_document:
    number_sections: true
    keep_tex: true
geometry: margin=0.5in
---

<!-- # Follow-up experiment (September 20, 2020) -->
# Preamble: Define functions, load data, and prepare data.

```{r}
# Preamble

rm(list=ls())

library(tidyverse)
library(vegan)
library(lme4)
library(gridExtra)
library(ggfortify)
library(car)
library(gllvm)

theme_set(theme_minimal())

res.bac = c("p", "o", "s", "a", "v")
all.bac = c(res.bac, "inv", "total.res", "total.all")


# Insert a ggplot layer under existing ones.
`-.gg` <- function(plot, layer) {
    if (missing(layer)) {
        stop("Cannot use `-.gg()` with a single argument. Did you accidentally put - on a new line?")
    }
    if (!is.ggplot(plot)) {
        stop('Need a plot on the left side')
    }
    plot$layers = c(layer, plot$layers)
    plot
}

# To format the scientific notation in the plots.
scientific_10 <- function(x) {
  ifelse(x==0, "0",
    parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x)))
  )
}

# Function to calculate derived metrics from the raw data.
prepare.data = function(Data, DetectionLimit=0){
  # Replace NA with 0 for calculating statistics and plotting data.
  Data[is.na(Data)] = 0
  # Replace "0"s with lower detection limit.
  Data[Data==0] = DetectionLimit
  # Transform time-bewteen-disturbance event to disturbance frequency.
  Data$disturbance = 1/Data$disturbance
  # Put in dilutions used in the experiments.
  Data = cbind(Data, plate_dilution = 10^-5,
    microcosm_volume_ml = 6,
    glycerol_dilution = 0.5,
    plated_volume_ul = 25) %>%
    mutate(total_dilution = plate_dilution * glycerol_dilution *
      plated_volume_ul / (microcosm_volume_ml * 1000))

  Data.prop.all = Data[,c(res.bac, "inv")]/rowSums(Data[,c(res.bac, "inv")])
  Data.prop.res = Data[,c(res.bac)]/rowSums(Data[,c(res.bac)])

  names(Data.prop.all) = paste(c(res.bac, "inv"), ".prop.all", sep="")
  names(Data.prop.res) = paste(c(res.bac), ".prop.res", sep="")

  # Combine together the different measures.
  Data = cbind(Data, Data.prop.all, Data.prop.res,
    simpson = diversity(Data[,res.bac], index="simpson"),
    shannon = diversity(Data[,res.bac], index="shannon"),
    alpha = rowSums(Data[, res.bac] > DetectionLimit),
    simpson.inv = diversity(Data[, c(res.bac, "inv")], index="simpson"),
    shannon.inv = diversity(Data[, c(res.bac, "inv")], index="shannon"),
    alpha.inv = rowSums(Data[, c(res.bac, "inv")] > DetectionLimit))
  Data = cbind(Data,
    # Hill number of order 1
    D1 = exp(Data[, "shannon"]),
    D1.inv = exp(Data[, "shannon.inv"]))
  Data = cbind(Data,
    # Scaled Hill number of order 1
    D1.scaled = (Data[, "D1"] - 1)/(length(res.bac)-1),
    D1.inv.scaled = (Data[, "D1.inv"] - 1)/(length(res.bac)+1-1))
  return(Data)
}


autoplot.custom = function(fit){
  autoplot(fit, alpha=0.5,
    smooth.colour="red",
    label.colour="red")
}

# ggplot settings.
pub.layer = list(
  theme(panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.title.y = element_text(size=12),
    axis.title.x = element_text(size=12),
    strip.background = element_blank(),
    strip.text = element_text(angle = 0, hjust = 0.0, size = 12),
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 9),
    legend.position= 'right'),
    labs(colour ='Invaded',
    fill ='Invaded',
    linetype ='Invaded',
    alpha ='Invaded',
    shape ='Invaded'))

# ggplot settings.
fig.layer = list(
  scale_color_manual(values=c("no"="black", "yes"="red")),
  scale_fill_manual(values=c("no"="black", "yes"="red")),
  scale_x_continuous(breaks=1/c(2,4,8,16), labels=c("2\nHigh", "4", "8", "16\nLow")),
  xlab("Pulse frequency (every x days)"),
  theme_bw(),
  guides(color=guide_legend(override.aes = list(alpha=1))))


# New data used to plot model predictions.
newdat = expand_grid(
  disturbance = seq(1/16, 1/2, length.out=100),
  invaded = c("no", "yes"),
  total_dilution = 1)

# Load the data.
data.diversity = read.csv("data.csv", comment.char="#") %>%
  prepare.data(Data=.)
```



# The effects of \inv{} and \freq{} on the effective number of species of the resident community.

Here we perform regressions on the scaled effective numbers of species of the resident community samples. Resident communities in our experiment can contain 1 to 5 microbial species. We measure and scale the effective numbers of species of the samples so the scaled numbers range from 0 (only 1 species) to 1 (equal abundance of all 5 species), using the following formula:
\begin{equation}
  \hat{D} = \frac{D-1}{D_{max}-1} = \frac{D-1}{4}
\end{equation}
where $D$ is the effective number of species and $D_{max}$ is the maximum possible number of species, i.e. 5.

We run logistic regressions with quasibinomial family on $\hat{D}$ to assess the effects of \inv{} and \freq{}:
  \begin{equation}
    \label{eq:FullDiv}
  logit(\hat{D}) \sim (\text{pulse} + \text{pulse}^2) \times \text{invasion},
  \end{equation}
where "pulse" is the pulse frequency (\textit{e.g.}, $1/16$ for the regime with a pulse every 16 d) and "invasion" is a categorical variable coded as 0 for non-invaded and 1 for invaded samples.

```{r}
fig.D1.scaled = ggplot(data=data.diversity,
  mapping=aes(x=disturbance, y=D1.scaled, color=invaded, shape=invaded, fill=invaded)) +
  geom_jitter(alpha=0.3, size=3.5, position=position_jitter(width=0.02, seed=1)) +
  geom_text(aes(label=alpha), position=position_jitter(width=0.02, seed=1), size=3) +
  scale_y_continuous(breaks=(1:5 - 1)/(4), labels=1:5) +
  fig.layer +
  pub.layer +
  ylab("Effective number of species")

# Logistic regression.
fit.D1.scaled = glm(D1.scaled ~ (disturbance + I(disturbance^2))*invaded,
  data=data.diversity, family="quasibinomial")

drop1(fit.D1.scaled, test="LRT")
# Quasibinomial, so no AIC output.
# The full model has the smallest deviance.

# Model diagnostics.
autoplot.custom(fit.D1.scaled)
# The Q-Q plot shows some deviance, but overall not too bad...


# Make a dataframe to calculate the model predictions +/- the standard errors.
pred = predict.glm(fit.D1.scaled, newdata=newdat, type="link", se.fit=T) %>%
  magrittr::extract(c("fit", "se.fit")) %>%
  as.data.frame() %>%
  cbind(newdat, .) %>%
  transform(p.fit=plogis(fit), p.lwr=plogis(fit-se.fit), p.upr=plogis(fit+se.fit))

# FIGURE 2
# Plot the model predictions onto the data points.
fig.D1.scaled -
  geom_line(data=pred, aes(y=p.fit, linetype=invaded), show.legend=F, size=0.75) -
  geom_ribbon(data=pred, aes(y=NULL, ymin=p.lwr, ymax=p.upr), show.legend=F, color=NA, alpha=0.2)
ggsave("D1scaled.png")


# The model output.
summary(fit.D1.scaled)

# Calculate McFadden's pseudo R-squared:
##  1 - (deviance)/null.deviance
with(fit.D1.scaled, 1 - (deviance)/null.deviance)
```

It looks like \inv{}, \freq{}, and \intx{} affect the resident diversity.



# The effects of \freq{} on invader fitness.

## Proportion of the invader.

Here we assess the effects of \freq{} on the proportion of the invader in the community.

```{r}
fig.inv.prop = ggplot(data=subset(data.diversity, invaded=="yes"),
 mapping=aes(x=disturbance, y=inv.prop.all)) +
 # geom_boxplot(aes(group=disturbance), outlier.shape = NA, lwd = 0.5,
 #   width = 0.03, show.legend = F, color="red") +
 geom_jitter(width=0.01, height=0, alpha=0.5, size=2.5, color="red", shape=17, fill="red") +
 ylab("Proportion of the invader") +
 ylim(c(0,1))+
 fig.layer +
 pub.layer


# First attempt at the logistic regression.
fit.inv.prop = glm(
  cbind(inv, (p+o+s+a+v)) ~ disturbance,
  data = data.diversity,
  subset = (invaded=="yes"),
  family = "binomial")

# Model diagnostics.
autoplot.custom(fit.inv.prop)
# The model fit doesn't look great. The residual plot shows possible quadratic effect.

# Plot the model.
# Make a dataframe to calculate the model predictions +/- the standard errors.
pred.inv.prop = predict.glm(fit.inv.prop, newdata=newdat, type="response", se.fit=T) %>%
  magrittr::extract(c("fit", "se.fit")) %>%
  as.data.frame() %>%
  cbind(newdat, .) %>%
  transform(upr=fit+se.fit, lwr=fit-se.fit)

fig.inv.prop -
  geom_line(data=pred.inv.prop, aes(y=fit),
    size=0.75, linetype=2, color="red") -
  geom_ribbon(data=pred.inv.prop, aes(y=NULL, ymin=lwr, ymax=upr),
    alpha=0.2, color=NA, fill="red")
# Issues with over- and under-estimating the data... need to include another term as expected.

# Second attempt at the logistic regression: include disturbance^2 as a covariate.
fit.inv.prop2 = glm(
  cbind(inv, (p+o+s+a+v)) ~ disturbance + I(disturbance^2),
  data = data.diversity,
  subset = (invaded=="yes"),
  family = "binomial")

# Model diagnostics
autoplot.custom(fit.inv.prop2)
# The Q-Q plot could be better.

drop1(fit.inv.prop2, test="LRT")
# This model is much better than the first one.

# Make a dataframe to calculate the model predictions +/- the standard errors.
pred.inv.prop2 = predict.glm(fit.inv.prop2, newdata=newdat, type="link", se.fit=T) %>%
  magrittr::extract(c("fit", "se.fit")) %>%
  as.data.frame() %>%
  cbind(newdat, .) %>%
  transform(p.fit = plogis(fit), p.upr=plogis(fit+se.fit), p.lwr=plogis(fit-se.fit))

# FIGURE 3A
# Plot the model predictions onto the data points.
fig.inv.prop -
  geom_line(data=pred.inv.prop2, aes(y=p.fit), size=0.75, linetype=2, color="red") -
  geom_ribbon(data=pred.inv.prop2,
    aes(y=NULL, ymin=p.lwr, ymax=p.upr), alpha=0.2, color=NA, fill="red")
ggsave("inv_prop.png")
# The fit looks much better!
# The model output.
summary(fit.inv.prop2)


# Calculate McFadden's pseudo R-squared:
##  1 - (deviance)/null.deviance
with(fit.inv.prop2, 1 - (deviance)/null.deviance)
```

It looks like \freq{} affects the proportion of the invader in the community.





## Density of the invader.

Here we assess the effect of \freq{} on the invader density.

```{r}
fig.inv.density = ggplot(data=subset(data.diversity, invaded=="yes"),
 mapping=aes(x=disturbance, y=inv/total_dilution)) +
 # geom_boxplot(aes(group=disturbance), outlier.shape = NA, lwd = 0.5,
 #   width = 0.03, show.legend = F, color="red") +
  geom_jitter(position=position_jitter(width=0.01, seed=21, height=0),
    alpha=0.5, size=2.5, color="red", shape=17, fill="red") +
  scale_y_continuous(label=scientific_10) +
  ylab("Invader density (CFU per microcosm)") +
  fig.layer +
  pub.layer +
  theme(legend.position="none")


# Attempt 1: Poisson regression.
fit.inv.density = glm(
  inv ~ disturbance + I(disturbance^2) + offset(log(total_dilution)),
  data=data.diversity,
  subset=(invaded=="yes"),
  family="poisson")

# Model diagnostics
autoplot.custom(fit.inv.density)
# Not too bad. The Q-Q plot could be better. Scale-Location plot shows some overdispersion?

drop1(fit.inv.density, test="LR")
# Better not to drop any term.

summary(fit.inv.density)
# Slightly overdispersed?


# Attempt 2: Negative binomial regression
fit.inv.density2 = MASS::glm.nb(
  inv ~ disturbance + I(disturbance^2) + offset(log(total_dilution)),
  data = data.diversity,
  subset = (invaded=="yes"),
  link = "log")

# Model diagnostics
autoplot.custom(fit.inv.density2)

AIC(fit.inv.density, fit.inv.density2)
# The negative binomial model is slightly better, but not by much...

# Add another term?
add1(fit.inv.density2, ~ . + I(disturbance^3), test="LRT")
# No.

# Plot the model.
pred.inv.density = predict.glm(fit.inv.density2, newdata=newdat, type="link", se.fit=T) %>%
  magrittr::extract(c("fit", "se.fit")) %>%
  as.data.frame() %>%
  cbind(newdat, .) %>%
  transform(upr=fit+se.fit, lwr=fit-se.fit) %>%
  transform(r.fit=exp(fit), r.upr=exp(upr), r.lwr=exp(lwr))

# FIGURE 3B
fig.inv.density -
  geom_line(data=pred.inv.density, aes(y=r.fit),
    size=0.75, linetype=2, color="red") -
  geom_ribbon(data=pred.inv.density, aes(y=NULL, ymin=r.lwr, ymax=r.upr),
    alpha=0.2, color=NA, fill="red")
ggsave("fig_density.png")

# The model output.
summary(fit.inv.density2)


# Calculate McFadden's pseudo R-squared:
##  1 - (deviance)/null.deviance
with(fit.inv.density2, 1 - (deviance)/null.deviance)
```

It looks like the \freq{} affects the density of the invader.









# Density of the community.

## The effects of \freq{} on the total (residents + invader) density

Here we test whether the \inv{} affects the total (resident + invader) density of the community.

```{r}
fig.total.density = ggplot(data=data.diversity,
    mapping = aes(x=disturbance, y = (p+o+s+a+v+inv) / total_dilution,
    color=invaded, shape=invaded)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.01,
      jitter.height = 0, dodge.width = 0.05, seed=1),
      alpha=0.5, size=2) +
  scale_y_continuous(label=scientific_10) +
  ylab("Resident + invader density (CFU per microcosm)") +
  fig.layer +
  pub.layer


fig.total.density
# There's a clear outlier, so remove it.
fig.total.density = fig.total.density %+%
  subset(data.diversity, (p + o + s + a + v + inv)/total_dilution <= 1e10)
# It looks like the total density changes with pulse frequency, but not with invasion.

# FIGURE 4B
fig.total.density -
  geom_boxplot(
    aes(group = interaction(disturbance, invaded)),
    outlier.shape = NA, lwd = 0.5,
    position = position_dodge(width = 0.05),
    width = 0.035, show.legend = F)
ggsave("total_density.png")


# Testing the effects of invasion as a two-way ANOVA.
fit.total.density = MASS::glm.nb(
  I(p+o+s+a+v+inv) ~ as.factor(disturbance) * invaded +
    offset(log(total_dilution)),
  data=data.diversity,
  subset=(p + o + s + a + v + inv)/total_dilution <= 1e10,
  link="log")

car::Anova(fit.total.density, test="LR")
# Invasion has slightly significant effect on total density (p=0.09)

# The model output.
summary(fit.total.density)
```

Invasion doesn't seem to affect the total community density.
Does an increase in invader density lead to a decrease in the resident density?




## The effects of \inv{} on the resident community density.

The analyses above suggest that the total densities do not differ between invaded and uninvaded communities. We will test this hypothesis below.

```{r}
fig.res.density = ggplot(data=data.diversity,
  mapping=aes(x=disturbance, y=(p+o+s+a+v)/total_dilution,
    color=invaded, shape=invaded)) +
  geom_boxplot(
    aes(group = interaction(disturbance, invaded)),
    outlier.shape = NA, lwd = 0.5,
    position = position_dodge(width = 0.05),
    width = 0.035, show.legend = F) +
  geom_point(position = position_jitterdodge(jitter.width = 0.01,
    jitter.height = 0, dodge.width = 0.05, seed=1),
    alpha=0.5, size=2) +
  ylab("Resident density (CFU per microcosm)") +
  scale_y_continuous(label=scientific_10) +
  fig.layer +
  pub.layer

fig.res.density
# One of the data points is clearly an outlier. Subset the data for the analyses.

fig.res.density = fig.res.density %+%
  subset(data.diversity, (p + o + s + a + v + inv)/total_dilution <= 1e10)
# The model that best captures the effects of disturbance on the community density is complex.
# So we will perform a simple two-way ANOVA to see if the invaded communities differ
# from uninvaded ones.

fit.res.density = MASS::glm.nb(
  I(p+o+s+a+v) ~ invaded*as.factor(disturbance) +
    offset(log(total_dilution)),
  data=data.diversity,
  subset=(p + o + s + a + v + inv)/total_dilution <= 1e10,
  link="log")

# Model diagnostics
autoplot.custom(fit.res.density)
# The Q-Q plot doesn't look too bad.

# Two-way ANOVA.
car::Anova(fit.res.density, test="LR")
# Significant effects of invasion and the interaction between invasion and disturbance.

# FIGURE 4A
fig.res.density
ggsave("res_density.png")
```

\Inv{} and \intx{} affect the resident density.










# Multivariate analyses on the individual resident species using Latent Variable Models.

Here we use the \texttt{gllvm} package to perform multivariate analyses on the density of each resident species.

```{r}
# library(gllvm)

# Labels for the species' names to be used in the plots.
labs = as_labeller(c(
    `p` = "P. corrugata",
    `o` = "O. daejonense",
    `s` = "S. rhizophila",
    `a` = "A. agilis",
    `v` = "V. guangxiensis"
))

# Create the data figure.
fig.indiv.res = data.diversity %>%
  subset((p + o + s + a + v + inv)/total_dilution <= 1e10) %>%
  pivot_longer(cols=p:v, names_to="spp", values_to="cfu") %>%
  ggplot(data=.,
    aes(x = disturbance, y = cfu/total_dilution,
      color = invaded, shape = invaded)) +
    geom_point(position = position_jitterdodge(jitter.width = 0.01,
      jitter.height = 0, dodge.width = 0.05, seed=1),
      alpha=0.3) +
    facet_wrap(facets = "spp", labeller = labs, nrow=1) +
    ylab("CFU per microcosm") +
    scale_y_continuous(label=scientific_10) +
    fig.layer +
    pub.layer +
    theme(strip.text = element_text(face = "italic", size=9),
      axis.text.x = element_text(size = 6.5))

data.temp = subset(data.diversity, (p + o + s + a + v + inv)/total_dilution <= 1e10)
data.temp$disturbance = as.factor(data.temp$disturbance)

# Analysis using factor(disturbance)

# FIGURE 5
# Boxplot
fig.indiv.res -
  geom_boxplot(
    aes(group = interaction(disturbance, invaded)),
    outlier.shape = NA, lwd = 0.5,
    position = position_dodge(width = 0.05),
    width = 0.035, show.legend = F)
    # Changing facet_wrap, otherwise hard to see.
    ggsave("indiv_density.png", width=12)
# Each resident species respond differently to invasion and pulse disturbances.

aicc.gllvm = rep(NA, 3)
gllvm.list = vector(mode="list", length=3)

for(i in 1:length(aicc.gllvm) ){
  fit = gllvm(
    y = dplyr::select(data.temp, p:v),
    X = data.temp,
    formula = ~ invaded * disturbance,
    offset = data.temp$total_dilution,
    num.lv = i-1,
    family = "negative.binomial", row.eff = "random",
    seed=1)

  gllvm.list[[i]] = fit
  aicc.gllvm[i] = summary(fit)$AICc
}

# Analysis with the model with the lowest AICc.
# Find the model with the lowest AICc.
which.min(aicc.gllvm)-1
# Number of latent variables: 0

# Set the model.
fit.gllvm0 = gllvm.list[[which.min(aicc.gllvm)]]

coefplot(fit.gllvm0, cex.ylab=1.2, cex.xlab=0.9,
  cex=1.25, order=F, mar=c(4, 3, 2, 1),
  mfrow=c(3,3))

summary(fit.gllvm0)
```






\appendix

\section*{Appendix: Additional analyses}

# The effects of \inv{} and \freq{} on the Gini-Simpson index of the resident community diversity.

The following figure shows the effects of \freq{} on the diversity of
the resident community, without (black circle) and with (red triangle) an \inv{}.
Light background points are individual sample measurements, and the dark points and their error bars are their mean and the standard error of the mean.
The number in each sample measurements is the species richness (i.e. the number of detected species; lower detection limit in this experiment: $`r 1/25/(10^-5)/0.5*1000`$ cfu per ml).

```{r}
fig.simpson = ggplot(data=data.diversity,
    mapping=aes(x=disturbance, y=simpson, color=invaded, shape=invaded)) +
  scale_color_manual(values=c("no"="black", "yes"="red")) +
  scale_fill_manual(values=c("no"="black", "yes"="red")) +
  scale_x_continuous(breaks=1/c(2,4,8,16), labels=c("2\nHigh", "4", "8", "16\nLow")) +
  geom_jitter(alpha=0.2, size=3.5,
    position=position_jitterdodge(dodge.width=0.05, jitter.width=0.01, seed=1)) +
  geom_text(aes(label=alpha), size=3, show.legend=F,
    position=position_jitterdodge(dodge.width=0.05, jitter.width=0.01, seed=1)) +
  xlab("Pulse frequency (every x days)") +
  ylab("Gini-Simpson diversity index") +
  guides(color=guide_legend(override.aes = list(alpha=1)))

fig.simpson -
  geom_boxplot(aes(group=interaction(disturbance, invaded)), outlier.shape=NA,
    position=position_dodge(width=0.05), width=0.02, show.legend=F)
```

It looks like the effect of disturbance on diversity is hump-shaped, and the effect of disturbance on diversity differs between invaded and uninvaded samples.

```{r}
lm.simpson1 = lm(simpson ~ (disturbance + I(disturbance^2)) * invaded, data=data.diversity)
summary(lm.simpson1);

# AIC comparison
drop1(lm.simpson1, test="F")

# Model with a lower AIC, but not significantly different from the full model.
lm.simpson2 = update(lm.simpson1, ~ . - I(disturbance^2):invaded)
summary(lm.simpson2)
```



```{r}
# Plot using the full model.
# Make a dataframe to calculate the model predictions +/- the standard errors.
pred = predict.lm(lm.simpson1, newdata=newdat, type="response", se.fit=T) %>%
  magrittr::extract(c("fit", "se.fit")) %>%
  as.data.frame() %>%
  cbind(newdat, .) %>%
  transform(lwr=fit-se.fit, upr=fit+se.fit)

# FIGURE S1A
fig.simpson -
  geom_line(data=pred, aes(y=fit, linetype=invaded), size=0.75) -
  geom_ribbon(data=pred, aes(y=NULL, ymin=lwr, ymax=upr, fill=invaded),
    show.legend=F, color=NA, alpha=0.1)+
  theme(legend.position="none")


# Plot using the model without I(disturbance^2):invaded.
# Make a dataframe to calculate the model predictions +/- the standard errors.
pred = predict.lm(lm.simpson2, newdata=newdat, type="response", se.fit=T) %>%
  magrittr::extract(c("fit", "se.fit")) %>%
  as.data.frame() %>%
  cbind(newdat, .) %>%
  transform(lwr=fit-se.fit, upr=fit+se.fit)

# FIGURE S1B
fig.simpson -
  geom_line(data=pred, aes(y=fit, linetype=invaded), size=0.75) -
  geom_ribbon(data=pred, aes(y=NULL, ymin=lwr, ymax=upr, fill=invaded),
    show.legend=F, color=NA, alpha=0.1)+
  theme(legend.position="none")
```