DataViz2.rmd

---
title: "Data Visualization 2"
author: "Joshua F. Wiley"
date: "`r Sys.Date()`"
output: 
  tufte::tufte_html: 
    toc: true
    number_sections: true
---

Download the raw `R` markdown code here
[https://jwiley.github.io/MonashHonoursStatistics/DataViz2.rmd](https://jwiley.github.io/MonashHonoursStatistics/DataViz2.rmd).


```{r loadpackages}

options(digits = 2)

## load relevant packages
library(tufte)
library(haven)
library(data.table)
library(JWileymisc)
library(ggplot2)
library(ggpubr)
library(ggthemes)
library(scales)
library(ggExtra)

## turn off some notes from R in the final HTML document
knitr::opts_chunk$set(message = FALSE)

```

```{r loaddata}

## read in data
db <- as.data.table(read_sav("B 19032020.sav")) # baseline
dd <- as.data.table(read_sav("DD 19032020.sav")) # daily


## average items then muliply to get back to "sum" scale
db[, PosAff := rowMeans(.SD, na.rm = TRUE) * 10,
   .SDcols = c("PANAS1", "PANAS3", "PANAS5", "PANAS9",
               "PANAS10", "PANAS12", "PANAS14", "PANAS16",
               "PANAS17", "PANAS19")]

```

# All Categorical Variables

If we are working with all categorical variables, a common way to
present them is to make one variable "continuous" by calculating the
percentages. For example, suppose that in our baseline data collection
exercise data, we wanted to graph the association between categorical
age and gender. 

```{r}

## create a categorical age variable
db[, AgeCat := factor(age < 22, levels = c(TRUE, FALSE),
                      labels = c("< 22y", ">= 22y"))]


egltable("AgeCat", g = "female", data = db)

```

## Bar Plot

After creating the categorical age variable, we can
get a frequency table easily enough (as shown in the earlier code),
but what if we wanted to graph it?

We could graph the frequencies, such as with a bar plot.

```{r}

p.bar <- ggplot(db, aes(female, fill = AgeCat)) +
  geom_bar(position = "dodge")

print(p.bar)

```

Beyond the data to ink ratio issues, if we were presenting this or
putting it into a paper, we would want to label it more cleanly.
Here we specify specific breaks on the x axis and specify their labels
and then remove the axis title, since saying men and women makes it
clear enough we don't need to say "female" the variable name anymore.
We could relabel the y axis (although count is fairly clear, I wanted
to show how to do it). The theme cleans up and makes the font a bit
bigger. 


```{r}

p.bar2 <- p.bar  +
  scale_x_continuous(
    breaks = c(0, 1),
    labels = c("Men", "Women")) +
  xlab("") +
  ylab("Frequency") +
  theme_pubr()

print(p.bar2)
 
```

`r margin_note("You can find more about customizing guides here: https://ggplot2.tidyverse.org/reference/guide_legend.html")`

We can change `AgeCat` by using 
`scale_fill_manual()` which lets us name the title of the legend and
to specify the colours, by name or hexademical codes, for each group
to make it black and white. We use the `coord_cartesian()` function to
stop the axis expansion so that it begins exactly at zero.
Finally, we get a title, with math symbols by using the `ggtitle()`
function listing the chi-square p-value from `egltable()` analysis
earlier. 

```{r}

p.bar3 <- p.bar2 +
  scale_fill_manual(
    "Age Group",
    values = c(
      "< 22y" = "black",
      ">= 22y" = "grey50")) +
  coord_cartesian(expand=FALSE) +
  ggtitle(expression(chi^2~p==.91))

print(p.bar3)

```

## Percentage Plot

A simple way would be to calculate the percentage of `female = 1` in
each age category and plot that. We would create a new dataset with
percentages calculated along with confidence intervals using the
following code. We calculate the average number of `female == 1` which
is the proportion, and then use the `prop.test()` function which takes
the count in one group and the total count and can calculate
confidence intervals to get 95% confidence intervals for the proportions.

```{r}

propdata <- db[!is.na(AgeCat), .(
  Percent = mean(female == 1, na.rm = TRUE),
  LL = prop.test(
    x = sum(female == 1, na.rm = TRUE),
    n = sum(!is.na(female)), correct = FALSE)$conf.int[1],
  UL = prop.test(
    x = sum(female == 1, na.rm = TRUE),
    n = sum(!is.na(female)), correct = FALSE)$conf.int[2]),
  by = AgeCat]

print(propdata)

```

Now we can plot the results. All we need for a basic plot is the
`ggplot()` and `geom_pointrange()` but the rest helps polish up the
figure for presentation.

```{r}

p.prop1 <- ggplot(propdata, aes(AgeCat, y = Percent, ymin = LL, ymax = UL)) +
  geom_pointrange() +
  scale_y_continuous(labels = percent) +
  scale_x_discrete(
    breaks = c("< 22y", ">= 22y"),
    labels = c("Age < 22y", "Age >= 22y")) + 
  xlab("") + ylab("Percent Female (95% CI)") + 
 theme_pubr()

print(p.prop1)

```

The same basic strategy can work for many variables at scale fairly
easily. For example, suppose that the personality dimensions were all
categorical. We first create these categorical variables, noting that
this is solely for the sake of demonstration. In general it is not a
good idea to convert continuous variables to categorical ones.

Next, we select just the variables we want (personality, ID, and age)
and reshape the dataset long where each variable is a
"timepoint". This allows us to have data table calculate the
proportions of each easily by setting by age category and by
personality variable. The resulting dataset has proportion and
confidence intervals of high on each personality measure for each age
group. 

```{r}

db[, O := as.integer(openness > median(openness, na.rm=TRUE))]
db[, C := as.integer(conscientiousness > median(conscientiousness, na.rm=TRUE))]
db[, E := as.integer(extraversion > median(extraversion, na.rm=TRUE))]
db[, A := as.integer(agreeableness > median(agreeableness, na.rm=TRUE))]
db[, N := as.integer(neuroticism > median(neuroticism, na.rm=TRUE))]

dblong <- reshape(
  db[!is.na(AgeCat), .(ID, AgeCat, O, C, E, A, N)],
  varying = list(Score = c("O", "C", "E", "A", "N")),
  v.names = "Score",
  timevar = "Personality",
  times = c("O", "C", "E", "A", "N"),
  idvar = "ID",
  direction = "long")

propdata2 <- dblong[, .(
  Percent = mean(Score == 1, na.rm = TRUE),
  LL = prop.test(
    x = sum(Score == 1, na.rm = TRUE),
    n = sum(!is.na(Score)), correct = FALSE)$conf.int[1],
  UL = prop.test(
    x = sum(Score == 1, na.rm = TRUE),
    n = sum(!is.na(Score)), correct = FALSE)$conf.int[2]),
  by = .(AgeCat, Personality)]

propdata2[, Personality := factor(Personality,
   levels = c("O", "C", "E", "A", "N"))]

print(propdata2)

``` 

Now we can plot the results. All we need for a basic plot is the
`ggplot()`, `geom_pointrange()`, and `facet_grid()` but the rest helps
polish up the figure for presentation. Note that we have not seen 
`facet_grid()` before. Facetting is an idea in data visualizing of
making "small multiples". Essentially the same plot over and over but
with some changes. In this case, its the same plot over and over but
changing the personality measure.

```{r}

p.prop2 <- ggplot(propdata2, aes(AgeCat, y = Percent, ymin = LL, ymax = UL)) +
  geom_pointrange() +
  scale_y_continuous(labels = percent) +
  scale_x_discrete(
    breaks = c("< 22y", ">= 22y"),
    labels = c("Age < 22y", "Age >= 22y")) + 
  xlab("") + ylab("Percent High (95% CI)") + 
  theme_pubr() +
  facet_grid(Personality ~ .) +
  coord_flip()

print(p.prop2)

```

# All Continuous Variables

For all continuous variables, there are not many graphing options. A
scatter plot is the main way two continuous variables are visualized.
However, even with a scatter plot, we can add additional
information to make it more useful. Our starting point is the scatter
plots with axes based on five number summaries we saw in Data
Visualization 1 topic. To that we add a linear regression line to help
show the overall association between the two variables. We also add a
text annotation with the correlation coefficient and p-value, found
from running the `cor.test()` function. We label it using `xlab()` and
`ylab()`. The main plot is saved in an object, `p.ss`.
Finally, we use the `ggMarginal()` function from the `ggExtra` package
to add histograms of the univariate distributions to the margins.
The final result captures extensive information about the individual
variables (through the histograms and five number summaries in the
axes) and about their association (through the scatter plot,
regression line, and correlation coefficient).


```{r, fig.height = 5.5, fig.width = 5.5, fig.cap = "scatter plot with regression line and correlation"}

cor.test(~ selfesteem + stress, data = db)

p.ss <- ggplot(db, aes(stress, selfesteem)) +
  geom_point() +
  stat_smooth(method = "lm", se = FALSE, size = 1) + 
  scale_x_continuous(breaks = as.numeric(quantile(db$stress))) + 
  scale_y_continuous(breaks = as.numeric(quantile(db$selfesteem))) +   
  theme_pubr() +   
  theme(axis.line = element_blank()) +
  geom_rangeframe() +
  xlab("Perceived Stress") +
  ylab("Self Esteem") +
  annotate("text", x = max(db$stress), y = max(db$selfesteem),
           label = "r = -0.54, p < .001",
           size = 6, hjust = 1, vjust = 1)

## now add a histogram to the margins
ggMarginal(p.ss, type = "histogram")

``` 

`r margin_note("Because we did not store the results of stress
and selfesteem with the histograms added to the margins, we have the
basic scatter plots in the figure, but we could have the histograms as
well if desired by saving the result in p.ss.")`

Another feature that is helpful with figures is to arrange sets of
related figures together. For example, in the following code we make a
plot of stress and neuroticism scores and save it in `p.sn`.
Now we can make a panel of graphs with two columns using the
`ggarrange()` function. 

```{r, fig.height = 5, fig.width = 10, fig.cap = "panel graph of two scatter plots"}

cor.test(~ neuroticism + stress, data = db)

p.sn <- ggplot(db, aes(stress, neuroticism)) +
  geom_point() +
  stat_smooth(method = "lm", se = FALSE, size = 1) + 
  scale_x_continuous(breaks = as.numeric(quantile(db$stress))) + 
  scale_y_continuous(breaks = as.numeric(quantile(db$neuroticism))) +   
  theme_pubr() +   
  theme(axis.line = element_blank()) +
  geom_rangeframe() +
  xlab("Perceived Stress") +
  ylab("Neuroticism") +
  annotate("text", x = max(db$stress), y = min(db$neuroticism),
           label = "r = 0.39, p = .003",
           size = 6, hjust = 1, vjust = 0)

ggarrange(
  p.ss, p.sn,
  ncol = 2,
  labels = c("A", "B"))

``` 

If we had more than two continuous variables we wanted to plot, there
are fewer options. We could create a 3D graph, but those rarely work
well for publications or theses that are printed and only viewed in
2D. Instead, its common to map additional variables to other aspects
or aesthetics in the figure. For example we can take our scatter plot
of stress and self esteem and make the size and shading of points
proportional to neuroticism scores. The following code does this and
plots the result. The only additions are adding
`size = neuroticism` and `colour = neuroticism` to 
`geom_point()`.

```{r, fig.height = 5.5, fig.width = 5.5, fig.cap = "scatter plot with regression line and correlation of stress, self esteem and neuroticism"}

cor.test(~ selfesteem + stress, data = db)

p.ssn <- ggplot(db, aes(stress, selfesteem)) +
  geom_point(aes(size = neuroticism, colour = neuroticism)) +
  stat_smooth(method = "lm", se = FALSE, size = 1) + 
  scale_x_continuous(breaks = as.numeric(quantile(db$stress))) + 
  scale_y_continuous(breaks = as.numeric(quantile(db$selfesteem))) +   
  theme_pubr() +   
  theme(axis.line = element_blank()) +  
  geom_rangeframe() +
  xlab("Perceived Stress") +
  ylab("Self Esteem") +
  annotate("text", x = max(db$stress), y = max(db$selfesteem),
           label = "r = -0.54, p < .001",
           size = 6, hjust = 1, vjust = 1)

print(p.ssn)

``` 

By default there are separate legend guides, one for size and one for
the colour. We can customize this using the `guides()` function in
`ggplot2`. We could turn one off or by using the same title make them
the same. Where we wrote `Neuroticism` capitalized, we also could have
changed the title of the legend guide (anything within the quotes
would be valid).

```{r, fig.height = 5.5, fig.width = 5.5, fig.cap = "scatter plot with regression line and correlation of stress, self esteem and neuroticism with one guide"}

p.ssn + 
  guides(
    size = guide_legend("Neuroticism"),
    colour = guide_legend("Neuroticism"))

``` 

# Mixed Continuous and Categorical Variables

The most variety of graphs are possible with data that are a mix of
continuous and categorical variables. As an example we will work with
age group, female, and the five personality traits. We start by
reshaping them long so that which personality measure is being
examined becomes another categorical variable.

```{r}

dblong2 <- reshape(
  db[!is.na(AgeCat), .(ID, AgeCat, female, openness, conscientiousness, extraversion, agreeableness, neuroticism)],
  varying = list(Score = c("openness", "conscientiousness", "extraversion", "agreeableness", "neuroticism")),
  v.names = "Score",
  timevar = "Personality",
  times = c("O", "C", "E", "A", "N"),
  idvar = "ID",
  direction = "long")
dblong2[, Personality := factor(Personality, levels = c("O", "C", "E", "A", "N"))]

head(dblong2)

``` 

We can make a simple plot with the means and 95% confidence intervals
using the following code.

```{r}

p.mean1 <- ggplot(dblong2, aes(Personality, Score)) +
  stat_summary(fun.data = mean_cl_normal) +
  theme_pubr()

print(p.mean1)

``` 

There are lots of additions we could consider to this simple
figure. The labels are short, but not the most informative. 
We could re-label the axis.

```{r}

p.mean2 <- p.mean1 +
  scale_x_discrete("",
    breaks = c("O", "C", "E", "A", "N"),
    labels = c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism"))

print(p.mean2)

```

Long labels are messy in smaller spaces. We could rotate the labels to
make space or rotate the graph.

```{r}

ggarrange(
  p.mean2 +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
    ggtitle("rotate text"),
  p.mean2 +
    coord_flip() +
    ggtitle("rotate graph"),
  ncol = 1)

```

Either rotating the labels or the graph made it easier to have long
labels and clearly read them, in this case especially rotating the
graph so the longest words can be read in their usual left to right
orientation.

With smaller datasets, we could show raw data as well as the means.
Although relatively easy to add, the result is disastrous. Even though
this is not a huge dataset, there is enough data and because there are
not many different possible values for each personality measure, the
dotplot is very difficult to read.
Lastly, because of the black dots and the means being shown as black
dots, it becomes impossible to well see the means.

```{r}

p.mean2 +
  geom_dotplot(binaxis = "y", stackdir = "center", binwidth = .2) + 
    coord_flip() +
    ggtitle("rotate graph")

```
By shrinking the size of the dots further (using `binwidth = .15`),
adding some transparency using the `alpha = .2` argument (valid
numbers are 0 completely transparent to 1 completely opaque) and
adding some random noise on the scores using `jitter()` we can see the
raw data and means better. It is not quite the actual raw data,
because we have added some noise, but it could still help to show the
general spread of data. 

```{r}

p.mean2 +
  geom_dotplot(aes(y = jitter(Score, 2)), binaxis = "y",
               alpha = .2,
               stackdir = "center", binwidth = .15) + 
    coord_flip()

```

For larger datasets or if the dotplot with noise is not as useful
anymore because its not true raw data, more summarized version of the
distribution can be shown with a violin plot, which is basically a
density plot that is mirrored. Thicker regions have more points,
narrow regions have fewer data points. We can also see the
range/spread of each variable and still have our mean and confidence
interval summaries shown clearly.

```{r}

p.mean3 <- p.mean2 +
  geom_violin(fill = NA) + 
  coord_flip()

print(p.mean3)

```

Another aspect we could imporve: there is no necessary ordering of
personality measures. Ordering them such as from highest to lowest
mean can be used to help us read the plot more easily.

We do this by changing the levels of the factor in the dataset.
First, we calculate the mean score by personality, then we have data
table order the resulting means from highest to lowest (these are
things we saw in the working with data topic).
Finally, we use `factor()` on personality and specify the levels in
this same order.

```{r}

dblong2[, .(M = mean(Score, na.rm = TRUE)), by = Personality][
  order(-M)]

dblong2[, Personality := factor(Personality,
                                levels = c("C", "A", "N", "O", "E"))]

```

Now we can simply remake our graph (note that this only works when
using data.table for data management, if using data frames etc. you
would want to copy and paste all your graph code again). 

By having ordered the variables by their means, it helps us rapidly
interpret which one has the lowest score (extraversion) and which the
highest average score (conscientiousness). It is a small step but one
that aids rapid processing of the figure and the data therein.

```{r}

print(p.mean3)

```

It is easy to add additional categorical variables into a figure.
For example, we can colour by age group.

```{r}

ggplot(dblong2, aes(Personality, Score, colour = AgeCat)) +
  stat_summary(fun.data = mean_cl_normal, position = position_dodge(.2)) +
  scale_x_discrete("",
    breaks = c("O", "C", "E", "A", "N"),
    labels = c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")) + 
  theme_pubr() +
  coord_flip()

```

If we wanted we could also add shapes by `female`. 
In this example, I also change the default colour and legend title for
`AgeCat` to Age Groups.

```{r}
    
ggplot(dblong2, aes(Personality, Score, colour = AgeCat, shape = factor(female))) +
  stat_summary(fun.data = mean_cl_normal, position = position_dodge(.3)) +
  scale_x_discrete("",
    breaks = c("O", "C", "E", "A", "N"),
    labels = c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")) + 
  theme_pubr() +
  coord_flip() +
  scale_colour_manual(
    "Age Group",
    values = c("< 22y" = "black", ">= 22y" = "grey80"))

``` 

If having four means side by side is too hard to read, we could facet
the plot into small multiples, say by female, so that we can compare
age groups in men and women.

```{r, fig.width = 7, fig.height = 10}

ggplot(dblong2, aes(Personality, Score, colour = AgeCat)) +
  stat_summary(fun.data = mean_cl_normal, position = position_dodge(.3)) +
  scale_x_discrete("",
    breaks = c("O", "C", "E", "A", "N"),
    labels = c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")) + 
  theme_pubr() +
  coord_flip() + 
  facet_grid(female ~ .) +
  scale_colour_manual(
    "Age Group",
    values = c("< 22y" = "black", ">= 22y" = "grey80"))

```

If we wanted to show the raw data, we could facet on both age and
female and add our dotplots back in.

```{r, fig.width = 10, fig.height = 9}

ggplot(dblong2, aes(Personality, Score)) +
  geom_dotplot(aes(y = jitter(Score, 2)), binaxis = "y",
               alpha = .2,
               stackdir = "center", binwidth = .35) +   
  stat_summary(fun.data = mean_cl_normal, position = position_dodge(.2)) +
  scale_x_discrete("",
    breaks = c("O", "C", "E", "A", "N"),
    labels = c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")) + 
  theme_pubr() +
  coord_flip() + 
  facet_grid(female ~ AgeCat)

```

When we facet, only the labels show up making it difficult to
interpret if this was being presented in presentation or article. In
this case, we might create new variables with more descriptive labels,
just for the plotting. Note that the choice of jitter, alpha, and
binwidth all involve some trial and error to get to a plot that is
easy to read and visually appealing (admittedly, a rather subjective
concept). 

```{r, fig.width = 10, fig.height = 9}

dblong2 <- copy(dblong2)
dblong2[, Sex := factor(female, levels = c(0, 1), labels = c("Men", "Women"))]
dblong2[, AgeCat2 := factor(AgeCat, levels = c("< 22y", ">= 22y"),
                            labels = c("Age < 22y", "Age >= 22y"))]

ggplot(dblong2, aes(Personality, Score)) +
  geom_dotplot(aes(y = jitter(Score, 2)), binaxis = "y",
               alpha = .2,
               stackdir = "center", binwidth = .35) +   
  stat_summary(fun.data = mean_cl_normal, position = position_dodge(.2)) +
  scale_x_discrete("",
    breaks = c("O", "C", "E", "A", "N"),
    labels = c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")) + 
  theme_pubr() +
  coord_flip() + 
  facet_grid(Sex ~ AgeCat2)

```

Next, we are going to look at some hypothetical data from an
intervention comparing augmented Treatment as Usual (TAU+) to
Cognitive Behavioural Therapy (CBT+). The two conditions are measured
at baseline and post intervention on depression symptoms.
The first part of the code just simulates some data including a wide
dataset, `trial` and a long dataset, `trial2`.
**You do not need to follow this code, it is just to get us some
sample data to work with.**

```{r}

## code to make an example dataset
set.seed(1234)
trial <- data.table(
  ID = sample(1:70),
  Group = factor(rep(c("TAU+", "CBT+"), each = 35)),
  B_Dep = pmax(round(rnorm(35*2, mean = 22, sd = 7)), 0))
trial[, P_Dep := round(B_Dep * rnorm(70, mean = ifelse(Group == "CBT+", .5, .9), sd = .2))]
trial2 <- reshape(trial, varying = list(c("B_Dep", "P_Dep")), v.names = "Depression",
                 timevar = "Assessment", times = c(0, 1),
                 idvar = "ID", direction = "long")

head(trial)

head(trial2)

```

With some sample data, we can plot the long dataset to show the mean
and confidence intervals for each group at each time point.

```{r}

p.trial1 <- ggplot(trial2, aes(Assessment, Depression, colour = Group)) +
  stat_summary(fun.data = mean_cl_normal,
               position = position_dodge(.05),
               geom = "pointrange") +
  theme_pubr()
p.trial1 <- set_palette(p.trial1, palette = "jco")

print(p.trial1)

```

Because these are longitudinal data, it makes sense to connect them
with lines to show how they changed over time. We do this by adding a
line geom based on the mean. Then we tidy up the x axis labels and the
y axis labels.

Finally, something new, we use `geom_hline()` to add a horizontal line
at 16, a common cut off on the CES-D indicative of clinically
significant depression symptoms. We make this a dashed, grey line to
make it less prominent. This line aids interpretation by helping
people anchor the results to common cut offs. We also use the 
`coord_cartesian()` function to change the limits of the graph. Since
the CES-D scale starts at 0 (meaning lowest possilbe / no depression
symptoms) we make that the y axis limit. The x axis limits are based
on the coding of assessments and the upper y axis limit we base
visually off the upper confidence interval.

```{r}

p.trial1b <- p.trial1 + 
  stat_summary(fun = mean,
               position = position_dodge(.05),
               geom = "line") +
  scale_x_continuous("",
                     breaks = c(0, 1),
                     labels = c("Baseline", "Post")) +
  scale_y_continuous("Depression Symptoms (CES-D)",
                     breaks = c(0, 4, 8, 12, 16, 20, 24)) + 
  geom_hline(yintercept = 16, linetype = 2, colour = "grey50") +
  coord_cartesian(xlim = c(-.05, 1.05), ylim = c(0, 26.5), expand = FALSE)

print(p.trial1b)

```

The other information that would be useful would be to annotate with
information about group differences and change over time.
First we run a regression on depression by group at each time point
and then use those p-values to add annotations to the graph.

```{r}

summary(lm(Depression ~ Group,
        data = trial2[Assessment == 0]))

summary(lm(Depression ~ Group,
        data = trial2[Assessment == 1]))

p.trial1b +
  annotate("text", x = 0, y = 26, label = "italic(n.s.)", parse = TRUE) + 
  annotate("text", x = 1, y = 26, label = "***")

```

In smaller datasets we could visualize the individual changes in
depression symptoms. We again plot depression symptoms on the y axis,
assessment on the x axis and colour by group, but instead of
summarizing the data, we directly plot points and lines. We use the
`group = ID` to indicate we want a different line for each ID in the dataset.

```{r}

p.trial2 <- ggplot(trial2, aes(Assessment, Depression, colour = Group, group = ID)) +
  geom_line() +
  geom_point() + 
  scale_x_continuous("",
                     breaks = c(0, 1),
                     labels = c("Baseline", "Post")) +
  scale_y_continuous("Depression Symptoms (CES-D)") + 
  geom_hline(yintercept = 16, linetype = 2, colour = "grey50") +
  theme_pubr()
p.trial2 <- set_palette(p.trial2, "jco")

print(p.trial2)

```

The result lets us see the starting point and change over time for
each person, but its a bit messy. Rather than just colour by group, it
might be helpful to separate by group, which we do by facetting.

```{r, fig.width = 9, fig.height = 6}

p.trial2 + facet_grid(. ~ Group)

```

That worked, but now our labels overlap. We need to add some space
between each facet (panel). Since each panel is labelled, we do not
really need the legend guide for group, so we turn that off by using
the `guides()` function to clean the plot up a bit.

```{r, fig.width = 9, fig.height = 6}

p.trial2 + facet_grid(. ~ Group) +
  theme(panel.spacing = unit(2, "lines")) +
  guides(colour = "none")

```

Another way to show individual change would be to use the wide dataset
to calculate individual change scores. A common approach is to examine
the percent change. We subtract 1 so that 0 means no change.
To plot the results, we put the individual IDs on the x axis and the
height of the bars is the percent change.

```{r}

trial[, PercentChange := P_Dep/B_Dep - 1]

p.trial3 <- ggplot(trial, aes(ID, PercentChange, fill = Group)) +
  geom_bar(stat = "identity") +
  theme_pubr() +
  scale_y_continuous("Change from Baseline", labels = percent)

p.trial3 <- set_palette(p.trial3, "jco")

print(p.trial3)

``` 

Although this figure is technically accurate, it is difficult to
interpret. The general pattern seems to be that the CBT+ group has a
more negative change. Ordering the data can improve this.
We use the `order()` function to order by percent change and then
order that to get numbers for a "new" ID variable. Now we can remake
the plot, with a few other tweaks to clean it up (a line at 0, no
change, removing the x axis and adding a better x axis title).

```{r}

trial[, ID2 := order(order(PercentChange))]

p.trial4 <- ggplot(trial, aes(ID2, PercentChange, fill = Group)) +
  geom_hline(yintercept = 0) + 
  geom_bar(stat = "identity") +
  theme_pubr() +
  scale_y_continuous("Change from Baseline", labels = percent) +
  xlab("Individual Participants") + 
  theme(
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank())

p.trial4 <- set_palette(p.trial4, "jco")

print(p.trial4)

```

With the data ordered, it is much easier to see the biggest decline,
the biggest increase and to see that the CBT+ group dominates the left
hand side with the largest decreases while little decrease or even
increases occur almost exclusively in the TAU+ group.

If we wanted a slight modification is to order first by group and then
by percent change, giving the followiing result.

```{r}

trial[, ID3 := order(order(Group, PercentChange))]

p.trial5 <- ggplot(trial, aes(ID3, PercentChange, fill = Group)) +
  geom_hline(yintercept = 0) + 
  geom_bar(stat = "identity") +
  theme_pubr() +
  scale_y_continuous("Change from Baseline", labels = percent) +
  xlab("Individual Participants") + 
  theme(
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank())

p.trial5 <- set_palette(p.trial5, "jco")

print(p.trial5)

```

# Summary Table

Here is a little summary of some of the functions used in this
topic. You might also enjoy this "cheatsheet" for `ggplot2`:
https://github.com/rstudio/cheatsheets/raw/master/data-visualization-2.1.pdf


| Function       | What it does                                 |
|----------------|----------------------------------------------|
| `ggplot()`     | Sets the dataset and which variables map to which aesthetics for a plot |
| `geom_point()` | Adds points such as for a scatter plot|
| `geom_hline()` | Adds a horizontal line at a specific y axis value |
| `stat_summary()` | Used to automatically calculate some summary statistics on data and plot, usually means with standard errors or confidence intervals | 
| `stat_smooth()` | Used to automatically calculate a regression line | 
| `ylab()` | Adds a label for the y axis |
| `xlab()` | Adds a label for the x axis |
| `theme_pubr()` | A cleaner black and white theme for `ggplot2` |