# Intro Case Study {#PG}

```{exercise, name="R Syntax"}

Start building a script by entering the commands into the text editor and then executing them in the console (press `ctrl` + `enter`). Commands will appear after the `>` sign.
```


```{r eval = FALSE}
PlantGrowth
```

`PlantGrowth` is a built-in data-set of plant dry-weight. If you execute `PlantGrowth` in the console, you'll see the entire data-set. The group column is a categorical variable, describing if dry-weight was measured under control (`ctrl`) or one of two different treatment (`trt1` and `trt2`) conditions. This experimental set-up is probably familiar to you, even if you are not studying plant growth.

We already have access to the data set, in the form of a `data frame`. To make it easier to print to the screen, we're going to convert it to a _tibble_.

```{r}
PlantGrowth <- as_tibble(PlantGrowth)
```

In many situations in statistics, _categorical variables_ are called _factors_. Within a factor, each _category_, or _group_, is called a _level_. We can find out how many levels we have and what they are by using the `nlevels()` and `levels()` functions on the growth column:

```{r}
nlevels(PlantGrowth$group)
levels(PlantGrowth$group)
```

Factors are very convenient because we can use them to split other columns. This is a very simple and powerful technique for handling data. To get a sense of the data's distribution, we can draw a variety of plots, sorting each weight value according to the levels of the growth factor (i.e. the three growth conditions). The result is displayed in figure \@ref(fig:boxplot2).

R specializes in handling and visualising data, so making the boxplot is straight-forward in `base` R. However, there are many new packages which have made things easier and more consistent, like `ggplot2`, which is included in the `tidyverse` suite. It offers fine control of making elegent statistical plots.

In the `ggplot()` function below, `weight` is `mapped` on the `x` scale, and `weight` is _mapped_ onto the `y` scale. Colloquially, we would way, "_weight_ against _group_", or "_weight_ as described by _group_":

```{r}
p <- ggplot(PlantGrowth, aes(group, weight))
```

Here is our box plot

```{r fig.cap = "A box plot", tidy = F}
p +
  geom_boxplot()
```

```{r tidy = F, echo = FALSE, eval = FALSE}
p +
  geom_jitter(width = 0.2, shape = 16, alpha = 0.75)
```

```{r tidy = F, echo = FALSE, eval = FALSE}
p +
  stat_summary(fun.data = mean_sdl, fun.args = list(mult = 1))
```

```{r fig.cap = "Summary statistics and raw values", tidy = F}
p +
  geom_jitter(width = 0.2, shape = 16, alpha = 0.75) +
  stat_summary(fun.data = mean_sdl, fun.args = list(mult = 1), col = "red")
```

We can get descriptive statistics on all the values in the `weight` column,

```{r}
mean(PlantGrowth$weight)
```

but group-wise statistics are more appropriate. You can imagine that statisticians have been doing such things for a long time. Let's use some functions that come with the `dplyr` package, also included in the `tidyverse` suite.^[You'll still see lots of old-school methods, namely the `apply` family of functions, from your colleagues and in online help forums.]

```{r tidy = F}
PlantGrowth %>%
  group_by(group) %>%
  summarise(avg = mean(weight),
            stdev = sd(weight))
```

As you can see from the example above, we can just keep adding to the `summarise` function -- that's really convenient! `%>%` is called the _pipe operator_, since we can _pipe_ many functions together, as shown here. When you read this function in your head, pronounce it as such: "Take `PlantGrowth` _and then_ group it by the `group` variable _and then_ summarize it into two columns: `avg`, which contains group-wise means and `stdev` which contains group-wise standard deviations of the `weight` values.  {#sec:pipe}

In the above plot, `stat_summary()` calculated and plotted the group-wise mean and standard deviation for us. Now that we've calculated those values ourselves, let's plot them directly. First, we need to save the output from above. If we just print it to the screen, there is no object to work on! Here, we'll use a little trick that works well in combination with the `%>%` operator: the "reverse assign" operator `->`. We only use this in the specific case of also using `%>%`.

```{r tidy = F, eval = TRUE}
PlantGrowth %>%
  group_by(group) %>%
  summarise(avg = mean(weight),
            stdev = sd(weight)) -> PlantGrowthSummary
```

Now I have an object that contains only my summary statistics. So let's plot it!

```{r fig.cap = "Summary statistics from their own, pre-calculated, data frame.", tidy = F}
ggplot(PlantGrowthSummary, aes(group, avg)) +
  geom_pointrange(aes(ymin = avg - stdev, ymax = avg + stdev))
```

You can probably guess what the new parts of this plotting function are doing. I didn't have to assign the _base_ layers of the plot to an object, `p` above, because I just used it once and then it was finished.

Notice that `PlantGrowthSummary` only has three rows, since calculating the mean is an aggregration function, i.e. it _aggergrates_ values to (typically) a single value. A Z-score is something quite different. This is a transformation, in which each value is _transformed_ in some way (i.e. _normalized_). Again `dplyr` offers a solution:

```{r tidy = F, eval = T}
PlantGrowth %>%
  group_by(group) %>%
  mutate(Z_score = scale(weight))
```

Applying a linear model is also straight-forward. Consider a one-way ANOVA:

```{r tidy = F}
Plant.lm <- lm(weight ~ group, data=PlantGrowth)
Plant.anova <- anova(Plant.lm)
```

This function calculate an ANOVA table on the linear model of `weight` described by `group` and saves it as an object called `Plant.anova`. If you type `Plant.anova` into the console, you will get a results table.

```{r comment="", prompt = TRUE, eval = TRUE, echo = FALSE}
Fvalue <- round(Plant.anova$`F value`[1],2)
p <- round(Plant.anova$`Pr(>F)`[1],3)
DF <- Plant.anova$Df
knitr::kable(Plant.anova, caption = 'The results of a one-way ANOVA on the plant growth data set.')
```

Since F (`r Fvalue`) is significant (p-value = `r p`) at the 5\% level on the F(`r DF`) distribution, we can reject the null hypothesis. This short workflow has already demonstrated R's three key strengths: data manipulation, statistical modelling, and data visualisation.

## Saving your Workspace

The R workspace refers to all objects you have defined in your current session. Functions for saving and loading the workspace, referred to as an _image_, are given in the following table.

Table: (\#tab:EnvFunc) Functions for managing your workspace.

| Area            | Function        | Description                                |
|-----------------|-----------------|--------------------------------------------|
| Workspace       | `ls()`          | List all objects in the workspace.         |
|                 | `rm(list=ls())` | Delete all objects in the workspace.       |
|                 | `save.image()`  | Save the whole workspace image to a file.  |
| Commands        | `history(n)`    | Display the last n entries in the history. |
|                 | `savehistory()` | Save history as a file.                    |
|                 | `loadhistory()` | Load a file as history.                    |
| Objects         | `saveRDS()`     | Save a single object to a file.            |
|                 | `save()`        | Save multiple objects to a file.           |
|                 | `load()`        | Load a workspace image or R object.        |
|                 | `dput()`        | Save an ASCII version of an R object.      |
|                 | `dget()`        | Load an ASCII version of an R object.      |

When saving your workspace, only objects in your environment are saved. The working directory, open graphics windows or loaded libraries (i.e. the current session state) are saved.

When you exit R, you will be asked to save the workspace image (i.e. environment) as `.RData`. **Don't do this!**. If you already exited RStudio and did happen to save the workspace image, the easiest thing to do is open the RStudio project but selecting it from the recent projects in the project pull down menu in the upper right corner of the IDE inteface or by opening (i.e. double-clicking) on the *.Rproj file in your project folder. If you have objects in your environment and haven't changed any of the default settings in RStudio, it came from the `.Rdata`. You may not see this file in your OS file listings unless you change some setting, since it's a hidden file (we can tell from the name, hidden files begin with a `.`) But it should be visible in the RStudio Files pane (lower-right). Select this file by clicking on the adjacent check-box, and delete it by clicking on Delete in the menu bar at the top of the Files pane.

The best practice is to start fresh with each R session and use well-commented scripts or markdown documents for managing your data analysis workflow. You may forget what you were working on by the next time you start R and, in the worst-case-scenario, you could have objects hanging around your environment which you're not even sure where they came from.

You'll probably encounter some of the following functions when looking for help online or if you inherit scripts from colleagues. They are outdated and should be avoided!

Table: (\#tab:AvoidFunc) Functions for managing your workspace that should be avoided.

| Area                 | Function        | Description                                |
|----------------------|-----------------|--------------------------------------------|
| Working Directory    | `getwd()`       | Check the working directory.               |
|                      | `setwd()`       | Change the working directory.              |
| Objects              | `attach()`      | "Attach" an object to the environment.     |
|                      | `detact()`      | "Detach" an object to the environment.     |