---
title: "08 An introduction to the package ggplot2"
output:
  html_document: default
  pdf_document: default
word_document: default
---
`emuR` has for historical reasons some specialized objects and specialized methods that allow working with these specialized `emuR` objects. While it is sometimes unavoidable to have such specialized objects and methods, it should be avoided to do so whenever possible - instead, we could use some standardized procedures that are very common in R.

In order to see the advantages of a more standardized procedure,
let us one again create a temporary EMU-SDMS-database first (using, by the way, some specialized (but unavoidable) commands from the package `emuR`):

```{r, echo=TRUE,results="hide",message=FALSE}
# load package
library(emuR)
# create demo data in directory
# provided by tempdir()
create_emuRdemoData(dir = tempdir())
# create path to demo database
path2ae = file.path(tempdir(), "emuR_demoData", "ae_emuDB")
# load database
ae = load_emuDB(path2ae, verbose = F)
```

#Specialized objects containing data need specialized commands to plot data

As we have seen in chapter 06, the default resulting object of a call to `get_trackdata()` is of class `trackdata`, which is a very special class only existing in the package emuR (and its predecessors).
The `emuR` package provides multiple specialized routines such as `dcut()`,
`trapply()`, `eplot` and `dplot()` for processing and visually inspect objects of this type (see Harrington, 2010, for the use of these functions). 
```{r}
vowels = query(ae,query="Phonetic==i:|u:|E")
vowels_fm = get_trackdata(ae,
                        seglist = vowels,
                        ssffTrackName =  "fm",
                        verbose = FALSE)
# show class of vowels_fm
class(vowels_fm)
```

The folloing command then extracts the formant values at the temporal midpoint of each segment (each vowel, in this case):
```{r}
vowels_fm05=dcut(vowels_fm,.5,prop = TRUE)
```

We can then use this object to plot the data and 95%-confidence ellipses.
```{r}
eplot(vowels_fm05[,1:2],label(vowels),centroid=TRUE,formant = TRUE)
```

The original `emutrack` trackdata object can be used to plot trajectories of formants (here: F2 only) as a function of time (first example) or a mean trajectory for each vowel categories' F2 as a function of normalized time (second example)

```{r}
dplot(vowels_fm[,2],label(vowels))
dplot(vowels_fm[,2],label(vowels),normalise=TRUE,average=TRUE)
```

These commands (and many other commands in the predecessors of `emuR`) are specialized to work with (and only with) `emutrack` data objects.

# Use of standard data.frames instead of specialized objects 

In most cases, however, a R user will store his data in data.frames. Data.frames are required for most commands in most packages concerned with plotting and/or statistical analyses.

As the `emutrack` trackdata object is a fairly complex nested matrix object with internal
reference matrices, which can be cumbersome to work with, the `emuR` package introduces
a new equivalent object type called emuRtrackdata that essentially is a
flat data.frame or data.table object. This object type can be retrieved by setting
the resultType parameter of the get trackdata() function to `emuRtrackdata`:

```{r}
vowels_fm_new = get_trackdata(ae,
                        seglist = vowels,
                        ssffTrackName =  "fm",
                        resultType="emuRtrackdata",
                        verbose = FALSE)
# show class of vowels_fm_new
class(vowels_fm_new)

vowels_fm_new

names(vowels_fm_new)
```

The emuRtrackdata object is an amalgamation of both a segment list and a trackdata object. The
first `sl_rowIdx` column of the `iVu` object indicates the row index of the
segment list the current row belongs to, the `times_rel` and `times_orig` (and `times_norm` in the forthcoming emuR-version) columns
represent the relative time and the original time of the samples contained in the
current row and `T1` (to `T`n in n dimensional trackdata) contains the actual signal
sample values. It is also worth noting that the `emuR` package provides a function
called `create emuRtrackdata()`, which allows users to create emuRtrackdata from
a segment list and a trackdata object. This is beneficial as it allows trackdata
objects to be processed using functions provided by the emuR package (e.g., dcut()
and trapply()) and then converts them into a standardized data.table object for
further processing (e.g., using R packages such as `lme4` or `ggplot2` which were implemented
to use with data.frame or data.table objects).

#Introduction to the package `ggplot2`

The goal of this chapter is to allow the reader to plot any numeric data from data.frames, whatever their source may be, including the new `emuRtrackdata` object. In order to do so, we sometimes have to manipulate the data.frame. We therefore will repeat some standard methods that manipulate data.frames.

The plots above can be done with ggplot2 and will look like:

![Figure 1: Equivalent to the eplot](Rgraphics/eplotrepl.png)

![Figure 2: Equivalent to the dplot](Rgraphics/dplotrepl.png)

![Figure 3: Equivalent to the normalized dplot](Rgraphics/dplotreplnorm.png)

Why `ggplot2`?
--------------

Advantages of ggplot2

-   consistent underlying `grammar of graphics` (Wilkinson, 2005)
-   plot specification at a high level of abstraction
-   very flexible
-   theme system for polishing plot appearance
-   mature and complete graphics system
-   many users, active mailing list

That said, there are some things you cannot (or should not) do With ggplot2:

-   3-dimensional graphics (see the rgl package)
-   Graph-theory type graphs (nodes/edges layout; see the igraph package)
-   Interactive graphics (see the ggvis package)

What Is The Grammar Of Graphics?
--------------------------------

The basic idea: independently specify plot building blocks and combine them to create just about any kind of graphical display you want. Building blocks of a graph include:

-   data
-   aesthetic mapping
-   geometric object
-   statistical transformations
-   scales
-   coordinate system
-   position adjustments
-   faceting

The structure of a `ggplot`
---------------------------

The `ggplot()` function is used to initialize the basic graph structure, then we add to it. The structure of a ggplot looks like this:

``` r
  ggplot(data = <default data set>, 
         aes(x = <default x axis variable>,
             y = <default y axis variable>,
             ... <other default aesthetic mappings>),
         ... <other plot defaults>) +

         geom_<geom type>(aes(size = <size variable for this geom>, 
                        ... <other aesthetic mappings>),
                    data = <data for this point geom>,
                    stat = <statistic string or function>,
                    position = <position string or function>,
                    color = <"fixed color specification">,
                    <other arguments, possibly passed to the _stat_ function) +

    scale_<aesthetic>_<type>(name = <"scale label">,
                       breaks = <where to put tick marks>,
                       labels = <labels for tick marks>,
                       ... <other options for the scale>) +

    theme(plot.background = element_rect(fill = "gray"),
          ... <other theme elements>)
```

The basic idea is that you specify different parts of the plot, and add them together using the `+` operator.

##Examples
See e.g. Handbook on R and figures: http://www.cookbook-r.com/
and the introduction to ggplot2 (gg = grammar of graphics) in http://docs.ggplot2.org/current/

Let's try with a few datasets from Jonathan Harrington's statistics seminar:

```{r}
# if necessary, install.packages(ggplot2)
library(ggplot2)

pfadu = "http://www.phonetik.uni-muenchen.de/~jmh/lehre/Rdf"
asp = read.table(file.path(pfadu, "asp.txt"))
coronal = read.table(file.path(pfadu, "coronal.txt"))
int.df = read.table(file.path(pfadu, "intdauer.txt"))
v.df = read.table(file.path(pfadu, "vdata.txt"))


# check class (data.frame or not):
class(asp)
# the first few lines:
head(coronal)


# 'ai[m,]' = row m
# 'ai[,m]' = column m
# You can use '$Name' to access column "Name"

#############################################################################
# 1. Numerical und categorical variables
############################################################################
# In a data.frame, columns can consist of numerical or categorical variables.
# In a matrix, you can only have one or the other class of variables.

# Numerical variables: continuous
#
class(asp$d)
# or
with(asp, class(d))
# [1] "numeric"

class(int.df$Dauer)
# [1] "integer"

# Categorical variables will be treated as factors (that have two or more levels, or categories; this is different to objects of the class "character"):
class(coronal$Socialclass)
# [1] "factor"

# first 10
coronal$Socialclass[1:10]
# asks which levels are given
levels(coronal$Socialclass)

##########################################################
# 2. Typical example in phonetics
##########################################################
# Is there an influence of x on y?
# 
# 1. y = numerical, x = categorical
# 1.1 difference in duration in /i, e, a/ ?
# 1.2 = influence of x (=vowel) on y (=duration)?
# 1.3 possible geoms: geom_boxplot() 
# or: geom_histogram()  or stat_density()

# 2. y = categorical, x = categorical
# 2.1 words like Sohn, Sonne... can be produced either with /s/ or /z/.
# /s/ more likely in Bavaria or in Hamburg?
# 2.2 possible geom: geom_barchart()

# 3. y = numerical, x = numerical 
# 3.1 bigger mouth opening related to a longer duration?
# 3.2 possible geom: geom_point(), geom_line()
```
### Summary: 
-   y = numerical, x = categorical: `geom_boxplot()`, `geom_histogram()`, `stat_density()`
-   y = categorical, x = categorical: `geom_bar()`
-   y = numerical, x = numerical: `geom_point()`, `geom_line()`
-   y = categorical, x = numerical: `geom_point()`

##Boxplots

```{r}
############################################################################
# 3. geom_boxplot(): y = numerical, x = categorical
############################################################################

head(asp)
# Influence of place of articulation (Kons) on duration of aspiration (d)?
# y: d    (numerical)
# x: Kons (categorical)

# Syntax in ggplot()
# A + B + C + D + ...
# A, B, C... are modules.
# Here:
# A. data-frame + B. Variables + C. kind of plot
ggplot(asp) + aes(y = d, x = Kons) + geom_boxplot()

# or
# A
p1 = ggplot(asp)
# B
p2 = aes(y = d, x = Kons)
# C
p3 = geom_boxplot()
# A + B + C
p1 + p2 + p3
# oder A + B + C ablegen
erg = p1 + p2 + p3
# Bild
erg

# boxplot.
# thick line = median;  'Box': interquartile range
# 
```

##Barplots

```{r}
############################################################################
# 4. geom_bar(): y ist kategorial, x ist kategorial
############################################################################
head(coronal)
# Influence of region (Region) in place of articulation (F1)?
# y: Fr (categorical)
# x: Region (categorical)

p1 = ggplot(coronal)
p2 = aes(fill = Fr, x = Region)
# to print frequencies of occurance
p3 = geom_bar()
p1 + p2 + p3
# place bars side by side
p4 = geom_bar(position="dodge")
p1 + p2 + p4

# print proportions
p5 = geom_bar(position="fill")
p1 + p2 + p5
```

##Scatterplots

```{r}
############################################################################
# 5. geom_point(), geom_line():  y ist numerisch, x ist numerisch
############################################################################  
# Inwiefern wird die Dauer (Dauer) von der Intensität (dB) beeinflusst in dem Data-Frame int.df()
# y: Dauer (numerisch) 
# x: dB (numerisch)
head(int.df)
# Nur Linie
ggplot(int.df) +  aes(x = dB, y = Dauer) + geom_line() 
# Nur Punkte
ggplot(int.df, aes(x = dB, y = Dauer)) + geom_point() 
# Beide
ggplot(int.df, aes(x = dB, y = Dauer)) + geom_line() + geom_point()
```

## Titles and axes

```{r}
############################################################################
# 6. + xlab() + ylab() + ggtitle()
############################################################################  
# same boxplot as above
p1 = ggplot(asp) + aes(y = d, x = Kons) + geom_boxplot()
# label for x-axis
p2 = xlab("Place of Articulation")
# label for x-axis
p3 = ylab("Duration (ms)")
# Titel
p4 = ggtitle("Boxplot")
p1 + p2 + p3 + p4

# same barchart as above
bar.p = ggplot(coronal) + aes(x = Region, fill = Fr) + geom_bar(position = "fill")
x.p = xlab("Region")
y.p = ylab("Proportion")
t.p = ggtitle("Proportional Distribution of Fricatives")
bar.p + x.p + y.p + t.p

############################################################################
# 7. Limits on axes +xlim() + ylim()
############################################################################

# same geom_bar() as above
p1 = ggplot(int.df, aes(dB, Dauer)) + geom_point() 
# xlim
p2 = xlim(c(10, 60))
# ylim
p3 = ylim(c(30, 280))
p1 + p2 + p3

#reverse axes:
p4 = scale_x_reverse()
p5 = scale_y_reverse()

p1 + p4 + p5
```




## Colors

(see  http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf)
```{r}
colors()

############################ geom_boxplot()
ggplot(asp) + aes(y = d, x = Kons) + geom_boxplot()
# Default colors
# filled with different colors
ggplot(asp) + aes(y = d, x = Kons, fill = Kons) + geom_boxplot()
# different line colors
ggplot(asp) + aes(y = d, x = Kons, col = Kons) + geom_boxplot()

# or chose your own colors
farben = c("green", "red")
# filled
ggplot(asp) + aes(y = d, x = Kons) + geom_boxplot(fill = farben)
# line colors
ggplot(asp) + aes(y = d, x = Kons) + geom_boxplot(col = farben)

############################ geom_bar()
########## 
p1 = ggplot(coronal) + aes(x = Region, fill = Fr) + geom_bar()
p1
# Eigene Farben wählen
farben = c("yellow", "green")
p2 = scale_fill_manual(values = farben) 
p1 + p2
```

## Plotting characters (pch) and character sizes (cex); line width (lwd)
(see http://www.endmemo.com/program/R/pchsymbols.php)

```{r}
########## 
ggplot(int.df, aes(x = dB, y = Dauer)) +  geom_point() + geom_line()
# col: color. 
# pch: plotting character. 
# cex: character expansion:cex =2 means 2*standard size
ggplot(int.df, aes(x = dB, y = Dauer)) +  geom_point(col="purple", pch=0, cex=2) + geom_line(col = "pink")
# lwd: Liniendichte
ggplot(int.df, aes(x = dB, y = Dauer)) +  geom_point(col="purple", pch=0, cex=2) + geom_line(col = "pink", lwd=2)
```



## Size of Labels

```{r}
# Default size ist 11 (Legende: 10 (??))

p1 = ggplot(asp) + aes(y = d, x = Kons) + geom_boxplot() + xlab("Artikulationsstelle") + ylab("Dauer (ms)") + ggtitle("Boxplot-Daten")
p1

# size 16
p16 = theme(text = element_text(size=16))
p1 + p16

# change only on axes
q24 = theme(axis.text = element_text(size=24))
p1 + q24

# Different values on axes labels and title
p30 = theme(text = element_text(size=30))
p1 + q24 + p30
```

## Two (or three) independent variables: Facets and (once again) colors

```{r}
#create one boxplot per stress pattern (Bet: levels "be" and "un")
pf = facet_grid(~Bet)
p1 + pf

# or add col to aes():
pc = ggplot(asp) + aes(y = d, x = Kons,col=Bet) + geom_boxplot() + xlab("Artikulationsstelle") + ylab("Dauer (ms)") + ggtitle("Boxplot-Daten")
pc
```

You can, of course combine facets and colors and therefore plot the influences of up to three independent variables.

## Arrange several figures on one plane

```{r}
# if necessary, install.packages(gridExtra)
library(gridExtra)

p1 = ggplot(asp, aes(y = d, x = Kons))  + geom_boxplot()
p2 = ggplot(coronal) + aes(x = Region, fill = Fr) + geom_bar()
p3 = ggplot(int.df, aes(dB, Dauer)) + geom_line() + geom_point()
grid.arrange(p1, p2, p3,  ncol=3, nrow =1)
```


## More things to change with `theme`
```{r, eval=FALSE}
# see
help(theme)
```

## Adding statistical measures to a plot

```{r}
p1 = ggplot(int.df, aes(dB, Dauer)) + geom_point()
int.lm =  geom_smooth(method="lm",se=FALSE)
p1 + int.lm
#by default, geom_smooth shows the standard error:
int.lmse =  geom_smooth(method="lm")
p1 + int.lmse

# you can calculate this stat (here lm() ) for each facet (e.g. for each subject (Vpn)) separately

p1 + int.lmse + facet_grid(~Vpn)


```

Instead of `geom_smooth()`, you could also add lines with `geom_abline(intercept=..., slope=... )`, and horizontal and vertical lines with `geom_hline()` and `geom_vline`.

`geom_smooth()` can be used with several smoothing methods, like `lm`, but also `glm` (for sigmoidal curves fitting binary perceptual data), and some others (it can fit e.g. splines with `loess`). One example of method `glm` (in which you have to add the information that it is binomial data) would be:

```{r, warning=FALSE}
bat.df = read.table("Rgraphics/dataSets/bat.df.txt")
bat.plot = ggplot(bat.df) + aes(y = p, x = steps) + geom_point(col = "red") +  facet_wrap(~participant)  + ggtitle("bat")

#add listener-specific sigmoids
bat.plot + geom_smooth(method = "glm",se=FALSE,method.args = list(family=binomial))

```

In phonetics, we often draw ellipses around two-dimensional data points, representing F2 and F1 values of vowels. We can add an ellipse by `stat_ellipse()`.

```{r}
ell = stat_ellipse()
p1 + ell
```

By default, this adds an ellipse representing the 95%-confidence interval (under the assumption of a multivariate t-distribution). While it is not extremely useful with the given data, it is useful in segregating vowel categories. However - be careful: at low numbers of tokens, one or two outliers can produce somehow "silly" ellipses: 

```{r}
td_mid = read.table("Rgraphics/dataSets/td_mid.txt")
p1 = ggplot(td_mid, aes(y = T1, x  = T2, col = labels, label=labels)) 

#add data.points as text labels, defined by their value
p2 = geom_text()
p1 + p2
p3 = stat_ellipse() 
p4 = scale_y_reverse() 
p5 = scale_x_reverse() 
p6 =labs(x = "F2(Hz)", y = "F1(Hz)") 
p7 = theme(legend.position="none")

p1 + p2 + p3 + p4 + p5 + p6 + p7

# only ellipses (do NOT plot data.points)
p1  + p3 + p4 + p5 + p6 + p7

#plot the label-specific means of F1 and F2 (here: T1 and T2)
p2_centroid = geom_text(data = aggregate(cbind(T1,T2)~labels,data=td_mid,FUN=mean))
p1 + p2_centroid + p3 + p4 + p5 + p6 + p7

#btw, we could also vary the linetype
p1_alt = ggplot(td_mid, aes(y = T1, x  = T2, col = labels, label=labels,linetype=labels)) 
p1_alt + p2_centroid + p3 + p4 + p5 + p6
```

It is also very easy to do the replacement of the dplot shown at the beginning of this document.
```{r}
ggplot(vowels_fm_new) +
  aes(x=times_rel,y=T2,col=labels,group=sl_rowIdx) +
  geom_line() +
  labs(x = "vowel duration (ms)", y = "F2 (Hz)")
```

However, it is much more difficult to produce the time-normalized and by-vowel averaged version.
We will need the function `normalizeLength()` (that will be available with the next release of `emuR`).
We can, however, use a prepared version of a length-normalized `emuRtrackdata` object that contains normalized times:

```{r}
td_norm = read.table("Rgraphics/dataSets/td_norm.txt")
ggplot(aggregate(T2~times_norm+labels, data = td_norm,FUN=mean)) +
  aes(x=times_norm,y=T2,col=labels) +
  geom_line() +
  labs(x = "vowel duration (normalized)", y = "F2 (Hz)")
```

## Conclusion
This chapter gave a very short introduction into the package ggplot2. More information can be found at e.g. http://r-statistics.co/Complete-Ggplot2-Tutorial-Part1-With-R-Code.html,
or https://opr.princeton.edu/workshops/Downloads/2015Jan_ggplot2Koffman.pdf,
or any other website you may find (there are numerous introductions to ggplot2).