
# Case Study

Effects of Meditation

To understand the usefulness of tidy data and `dplyr`, consider the _meditation_ data-set. This data-set contains the expression values of 6 pro-inflammatory genes. The researchers wanted to test if these genes are differentially expressed after meditation in people who regularly meditate, compared to non-meditating control individuals. Expression levels were measured in the morning (time 1) and after (time 2) a full day of meditation, by experienced meditators, or leisure, by the control group.


```{r excel, echo=FALSE, fig.cap="The `medi` data-set as it would be typically stored as an Excel spreadsheet. Although this is _human-readable_, it is _not_ an efficient way of storing data!", out.width = '50%'}
knitr::include_graphics("figures/Excel_Example_trans.jpg")
```

```{exercise, name="Messy versus Tidy"}
What is the problem with the data format in figure \@ref(fig:excel)?
```


Figure \@ref(fig:excel) is an example of a typical Excel spreadsheet. Here, column headers are found in the middle of the spread sheet and the results for each measured gene is contained in its own table on the same worksheet. In addition, column headers are colour-coded and there is plenty of white space for formatting purposes.


As a first step towards tidy data, we will remove all white space, colour-coding and floating labels. The data-set in Figure \@ref(fig:ExpBetter) resembles the long list file we have already encountered. It's not perfect, but it's a good start for working with R.

```{r ExpBetter, echo=FALSE, fig.cap="The `medi` data-set in a format ready for importing into R.", out.width = '100%'}
knitr::include_graphics("figures/ExpBetterFormat.jpg")
```

```{r prompt = FALSE, echo=FALSE, eval=TRUE, warning=FALSE, error=FALSE, message=FALSE}

read.delim("data/Expression.txt") %>% 
  tbl_df() -> medi

```

```{exercise, name="Import data"}
Import `Expression.txt`. Save it as an object called `medi` and convert it to a tibble.
```
Confirm that `medi` has the dimensions `r nrow(medi)` x `r ncol(medi)`.

<!-- We can aggregrate our \emph{values} (in this case take the mean) for each \emph{variable} according to \emph{id}, as per \ref{case4}, using the command: -->

<!-- without aggregration, case 4 -->

<!-- ```{r echo = TRUE, eval = FALSE} -->

<!-- dcast(df.m, type ~ variable, mean) -->

<!-- ``` -->


<!-- ```{r echo = FALSE, results = 'asis', eval = FALSE} -->

<!-- df.table <- xtable(dcast(df.m, type ~ variable, mean), caption = "The mean of each \\texttt{variable} category, separated by \\texttt{id}.", label = "case4") -->
<!-- digits(df.table) <- 0 -->
<!-- print(df.table, include.rownames=FALSE, floating.environment='margintable', latex.environments = "", table.placement = NULL) -->


<!-- ``` -->

<!-- Likewise, we can aggregrate our _values_ for each _variable_ according to _time_, as in \@ref(case5). -->

<!-- without aggregration, case 5 -->

<!-- ```{r echo = TRUE, eval = FALSE} -->

<!-- dcast(df.m, time ~ variable, mean) -->

<!-- ``` -->


<!-- ```{r echo = FALSE, results = 'asis', eval = FALSE} -->


<!-- df.table <- xtable(dcast(df.m, time ~ variable, mean), caption = "The mean of each \\texttt{variable} category, separated by \\texttt{time}.", label = "case5") -->
<!-- digits(df.table) <- 0 -->
<!-- print(df.table,include.rownames=FALSE, floating.environment='margintable', latex.environments = "", table.placement = NULL) -->

<!-- ``` -->


<!-- Case 6 -->

<!-- An alternative would be to aggregrate our \emph{values} for each \emph{time} according to \emph{id}, as in \ref{case6}. -->

<!-- without aggregration, case 6 -->

<!-- ```{r echo = TRUE, eval = FALSE} -->

<!-- dcast(df.m, type ~ time, mean) -->

<!-- ``` -->


<!-- ```{r echo = FALSE, results = 'asis', eval = FALSE} -->
<!-- df.table <- xtable(dcast(df.m, type ~ time, mean), caption = "The mean of each \\texttt{time} category, separated by \\texttt{id}.", label = "case6") -->
<!-- digits(df.table) <- 0 -->
<!-- print(df.table,include.rownames=FALSE, floating.environment='margintable', latex.environments = "", table.placement = NULL) -->
<!-- ``` -->

<!-- Let's return to the `medi` data-set. -->

If no ID variables are provided to `gather()`, it tries to guess which variables are the best choice. To do this, `gather()` uses categorical variables and then any column headers of continuous data as ID variables. However, recall that the `medi` data-set _doesn't have any_ categorical variables, so in this case `gather()` will take _all_ our column headers as a single ID variable.<div class="tooltip">&#9757;<span class="tooltiptext">The new column headings are what we defined in `gather()`: `key` and `value`.</span></div>

```{r}

# Original data-set:
dim(medi)

# Tidy data-set
medi.t <- gather(medi, key, value)
glimpse(medi.t)

```

The last action we need to perform before we can start working on this data is to separate the three variables contained in the `variable` column into distinct columns. For this we will use the `tidyr::separate()` function. Split up `medi.t$key` and overwrite `medi.t` with the new columns and the original values.

You should have a data frame with the following properties:<div class="tooltip">&#9757;<span class="tooltiptext">Note that the number of rows of `medi.t`, `r dim(medi.t)[1]`, is the product of the original data-set's dimensions, `r dim(medi)[1]` and `r dim(medi)[2]`.</span></div>



```{r echo = FALSE}
medi.t <- separate(medi.t, key, c("treatment", "gene", "time"), "_")
glimpse(medi.t)
```


<!-- % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EXERCISE Begin -->
<!-- % % \marginnote[1.5cm]{MARGIN} -->
<!-- % \marginnote[1.5cm]{\color{orange}\texttt{tidyr::spread()}} -->
<!-- % \begin{leftbar}{orange} -->
<!-- % \vspace{-1em} -->
<!-- % \subsection*{Exercise \theexercise} -->
<!-- % \stepcounter{exercise} -->
<!-- % \noindent{}Now that the \texttt{medi} data frame is tidy and each variable has its own column, we have access to all our variables. Generate the following table, which takes the mean of each \texttt{treatment} and \texttt{time} for every \texttt{gene}, using the \texttt{spread()} function: -->
<!-- % \bigskip -->
<!-- % \end{leftbar} -->
<!-- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EXERCISE End -->
<!-- %% <<echo = FALSE, results = 'asis'>>= -->
<!-- %% df.table <- xtable(dcast(medi.m, time + treatment ~ gene, mean, na.rm = TRUE)) -->
<!-- %% digits(df.table) <- 3 -->
<!-- %% print(df.table,include.rownames=FALSE) -->
<!-- %% @ -->
<!-- %% -->

## A Review of the 5 `dplyr` Verbs

Now that we have tidy data, let's review some of the functions from the \texttt{dplyr} package:

**filter** subsets the data frame. Alternatively, we could have used the base package function `subset()` or `[]`.

```{r}
medi.t2 <- filter(medi.t, value <= 1)
dim(medi.t2)

medi.t2 <- filter(medi.t, gene %in% c("CXCR1","RIPK2"))
dim(medi.t2)

medi.t2 <- filter(medi.t, !is.na(value))
dim(medi.t2)

```

**select** chooses specific variables

```{r}

medi.t2 <- select(medi.t, gene)
names(medi.t2)

medi.t2 <- select(medi.t, -gene)
names(medi.t2)

```


We can also rename variables as we select them:

```{r}
# Rename variables:
# Remember, select() keeps only the variables specified:
medi.t %>% 
  select(expressionValues = value)
```

```{r}
# rename() keeps all variables:
medi.t %>% 
  rename(expressionValues = value)
```

**arrange** reorders variables, similar to `order()`.

```{r}

# reorder according to descending value
medi.t %>% 
  arrange(desc(value)) -> medi.t2

# ascending order of gene, then time, then treatment
medi.t %>% 
  arrange(gene, time, treatment)

```

**mutate** creates new variables that are functions of existing variables

```{r}
# calculate the sd for the data set.
medi.t %>%
  mutate(stdev = sd(value, na.rm = T))
```

Alternatively, we can use `transmute()` which will drop existing variables:

```{r}
# calculate the sd for the data set.
medi.t %>%
  transmute(stdev = sd(value, na.rm = T))
```

**summarise** collapses multiple values to a single value, like in **ddply()**

```{r}
medi.t %>%
  summarise(average = mean(value, na.rm = T))
```

## The `group_by()` Adverb

So far so good, but of course in split-apply-combine, we are interested in working on sub-sets of our data frame. Therefore, in addition to the five nouns, there is one adverb, which modifies how the verb interacts with the noun. This is the important grouping variable:

`group_by()` defines the factor variables to use for subsequent actions.

```{r echo = FALSE}
options(tibble.print_max = 2)
```

```{r eval = FALSE}
# Treatment and time, grouped according to each gene:
medi.t %>%
  group_by(gene)

# Gene and time, grouped according to each treatment:
medi.t %>%
  group_by(treatment)

# Gene and treatment, grouped according to each time:
medi.t %>%
  group_by(time)

# Time, grouped according to each gene and treatment:
medi.t %>%
  group_by(gene, treatment)

# Treatment, grouped according to each gene and time:
medi.t %>%
  group_by(gene, time)

# Genes, grouped according to each treatment and time:
medi.t %>%
  group_by(treatment, time)

# All 24 combinations, grouped according to gene, treatment and time:
medi.t %>%
  group_by(gene, treatment, time)
```

The _grouped by_ object is still a data frame, but it has a new class attribute associated with it, which tells other `dplyr` functions how to group sub-sets.
```{r}
medi.t %>%
  group_by(gene, treatment, time) %>%
  class()
```

```{exercise, name="Calculate Statistics"}
Consider the structure of the data frame below. Using \texttt{dplyr} functions calculate each of the following statistics for each of the unique 24 combinations of gene, treatment and time:

- `average` the mean of the value
- `n` the number of observations in each group
- `SEM` The standard error of the mean
- `CIerror` The 95% CI error defined by the t distribution
- `lower95` The upper 95% CI limit
- `upper95` The upper 95% CI limit
```


```{r tidy = FALSE, echo = FALSE}

medi.t %>%
  group_by(gene, treatment, time) %>%
  filter(!is.na(value)) %>%
  summarise(average = mean(value),
            n = n(),
            SEM = sd(value)/sqrt(n),
            CIerror = qt(0.975, df = n-1),
            lower95 = average - (CIerror * SEM),
            upper95 = average + (CIerror * SEM)) -> medi.summary

library(DT)
datatable(medi.summary) %>% 
  formatRound(columns=c('average', 'SEM', 'CIerror', 'lower95', 'upper95'), digits=2)

```


```{r mediSumPlot, echo = FALSE, fig.cap = 'Gene expression averages and 95 CIs.', fig.width = large_width*1.5}
ggplot(medi.summary, aes(x = as.numeric(time), y = average, colour = treatment)) +
  geom_line(position = position_dodge(0.2)) +
  geom_pointrange(aes(ymin = lower95, ymax = upper95), position = position_dodge(0.2)) +
  facet_grid(. ~ gene) +
  scale_x_continuous(breaks = 1:2, labels = c("8:00", "16:00")) +
  labs(x = "Time", y = "Average Fold-Change gene expression (95CI)", col = "Treatment")
```

```{r mediRawPlot, echo = FALSE, fig.cap = 'The raw values.', fig.width = large_width*1.5}
# make a plot:
ggplot(medi.t, aes(time, value, col = treatment)) +
  geom_point(position = position_jitterdodge(0.2, dodge.width = 0.3)) +
  facet_grid(. ~ gene) +
  scale_x_discrete(breaks = 1:2, labels = c("8:00", "16:00")) +
  labs(x = "Time", y = "Fold Change", col = "Treatment")

```

```{exercise, name="Export data"}

Now that you\'ve processed your data, refer back to the chapter on importing data and see if you can get a file on your computer that contains the summary statistics. See section \@(exporting).

```
 
```{r eval = FALSE, echo = FALSE}
# save to your computer:
# use write.table(), write.csv(), or write.delim()
# or.... use the rio package
library(rio)
export(medi.summary, "summary data.txt")

```

```{exercise, name="Calculate ANCOVA"}

The `broom` package contains two functions to make applying statistical tests in a tidy format easy: `do()` and `tidy()`. See if you can calculate an ANCOVA for each gene in one command. Remember we saw how to calculate an ANOVA in the first chapter, an ANCOVA is a variant with a continuous predictor variable (here time). Out model will be defined as `value ~ treatment + time`.

This is the result for only RIPK2:

```

```{r echo = FALSE}
medi.t %>% 
  filter(gene == "RIPK2") %>% 
  lm(value ~ treatment + time, data = .) %>% 
  anova() %>% 
  kable()
```

In the end you should have a table that reports all the following values:

```{r ANCOVA, echo = FALSE}

# Do a single ANCOVA: use + or * for interaction
# use the broom package to clean up the results
library(broom)
# so... generalise to all genes: replace filter with group_by
medi.t %>% 
  # separate(key, c("treatment", "gene", "time"), "_") %>% 
  group_by(gene) %>% 
  do(tidy(anova(lm(value ~ treatment + time, data = .)))) %>% 
  datatable() %>% 
  formatRound(columns=c('sumsq', 'meansq', 'statistic', 'p.value'), digits=2)

```

```{exercise}
Produce a table reporting only the interesting values:
```

```{r ANCOVA2, echo = FALSE}
# so... generalise to all genes: replace filter with group_by
medi.t %>% 
  # separate(key, c("treatment", "gene", "time"), "_") %>% 
  group_by(gene) %>% 
  do(tidy(anova(lm(value ~ treatment + time, data = .)))) %>% 
  select(gene, term, p.value) %>% 
  filter(term != "Residuals") %>% 
  datatable() %>% 
  formatRound(columns=c('p.value'), digits=2)
```


```{exercise}
As a final exercise convert your script into an interactive markdown document so that you and your colleagues can explore the data tables more easily
```
