# Descriptive Statistics -- Understanding Your Data {#descstats}

Descriptive statistics are often given a cursory treatment, both in classes and in real--life analyses. People think of them as self--explanatory and somewhat boring and are impatient to jump to more sophisticated downstream methods, like hypothesis testing. In this chapter, we'll argue that using descriptive statistics is necessary for:

- Developing an intuitive feeling for your data,
- Finding mistakes in your data, 
- Generating hypotheses,
- Communicating results in presentations and publications (i.e. data visualization).

The reason why using descriptive statistics is so helpful with all of these things is that it provides a way to compress the dataset and to present it in a way that you and your colleagues can immediately understand. In essence there are two ways to carry out a descriptive analysis, either using one (or a few) numbers that describe the variables in an informative way or using plots. In this chapter, we will briefly talk about both of these options. The Scavetta Academy also offers an entire workshop dedicated to data visualization.

## Types of Variables {#variabletypes}

Before we start with the actual data description, one of the most important points of this chapter is to remind ourselves of the fact that _there are different kinds of variables_! To do so, let's take a look at the dataset in Table \@ref(tab:dataexcerpt) which shows some of the variables that were recorded for the 10 Martians collected at site I:

```{r dataexcerpt, echo = FALSE}
martian %>% 
  filter(Site == "Site I") %>% 
  select(-Group, -Site, -Height_UnequalVar8, -Height_UnequalVar, BMI, -Ratio, -AgeIndex, -Eye) %>% 
  knitr::kable(caption = 'Some variables measured on site I Martians.')
```



```{exercise, name="Different kinds of variables"}
What different kinds of variables do you see in Table \@ref(tab:dataexcerpt)?
```

## Categorical variables {#catdata}

Categorical variables describe ^[You may also see this referred to as _discrete_ data, which refers to distinct, non-overlapping groups. Other names are _qualitative_ or _factor_ data. Don't let the different terminology confuse you!] qualitative features and can take only a finite and countable set of possible values. In Table \@ref(tab:dataexcerpt), this description applies to the variables nose colour and eye sight. Despite this common feature, there is a clear difference between the two variables. There is no inherent order to the values of nose colour, they can't be put into a logical ascending or descending sequence. nose colour is a _nominal_ variable. In contrast, the values of eye sight can; they are ordered (poor < average < excellent). Eye sight is an _ordinal_ variable.


Table: (\#tab:cattypes) Types of categorical variables.

|         | Nominal                    | Ordinal                    |
|---------|----------------------------|----------------------------|
|         | ![](figures/Nominal4x.png) | ![](figures/Ordinal4x.png) |
| Ordered | no                         | yes                        |
| Example | nose colour                | eye sight                  |

## Continuous variables {#contdata}

Continuous variables describe quantitative features and can take any value on a continuous scale. Obvious examples in the Martian dataset are height and age. Like with categorical variables, there are two sub-categories that we distinguish for continuous variables, namely _interval_ and _ratio_.


Table: (\#tab:conttypes) Types of continuous scales.

|                 | Interval                    | Ratio                         |
|-----------------|-----------------------------|-------------------------------|
|                 | ![](figures/Interval4x.png) | ![](figures/Ratio4x.png)      |
| Meaningful zero | No                          | Yes                           |
| Example         | temperature in celsius      | height, age, temperature in K |

The difference between these sub--categories is rather academic: variables on a ratio scale have a natural zero whereas the value zero doesn't have a direct natural meaning for variables on an interval scale. To demonstrate this distinction, think about temperature. On the Kelvin scale, the value zero really refers to the absolute minimum, whereas $0$ is a rather arbitrary value in terms of absolute temperature. ^[The presence of a natural zero is important when it comes to building ratios. In terms of absolute temperature, 20K is really twice as warm as 10K whereas the relationship between 10 C and 20 C is not immediately obvious.]

As the vast majority of continuous variables will come with a natural zero, we won't bother with this distinction and will simply use the term _continuous_ for both interval and ratio variables for the rest of the workshop.

## Grey Zones

As always in life, things are not clear-cut, and sometimes it can be difficult to determine the variable type. Consider the number of antennae. From the dataset we see that Martians can have 0, 1, 5 or 9 antennae. Does this make the number of antennae continuous? Or is it ordinal with four categories? And now that we're at it, what if the eye sight would have been recorded in dioptre (a unit of refractive power)? Would that turn eye sight into a continuous variable?

Actually it is hard to think of an ordinal variable that isn't just a categorised version of a continuous variable. Whether you treat the variable as ordinal or continuous depends on the nature of the categorisation. With just three categories, eye sight is clearly an ordinal variable -- especially as we don't know how the categories were built and if there is an equal spacing between them. For number of antennae, the observed values in Table \@ref(tab:dataexcerpt) would also point towards an ordinal variable. But imagine the next Mars expedition results in a big dataset in which the number of antennae takes most values between 0 and e.g. 107, then treating number of antennae as continuous might be justified.

```{exercise, name="Ordinal variables"}
Can you think of situation where age would be treated as a ordinal variable?
```

## Suitable Questions and Information Content

The nature of the variables in a dataset is so important because it determines the questions that you can ask! This dictates the types of mathematical operations that can be performed which in turn has a direct impact on the type of analyses (including statistical tests) that can be used to investigate the data. Here is a summary of the types of questions that are meaningful for the different variable types that we covered above:


Table: (\#tab:questionstypes) The types of variable determines what questions we can ask!

|                                                          | Nominal | Ordinal | Interval | 
|----------------------------------------------------------|---------|---------|----------|
| How many of a certain kind?                              | yes     | yes     | yes      |
| How many values are bigger than a certain value?         | no      | yes     | yes      |
| What's the (mathematical) difference between two values? | no      | no      | yes      |

The above table tells us that information content increases from nominal to ordinal to continuous. The more information you have, the more questions you can ask of your data!

## Descriptive Summary Statistics

At the beginning of the chapter we said that the primary aim of descriptive statistics is to compress the dataset that you have in front of you in a way that the main features can be immediately understood. For a single variable, this is usually carried out by calculating two important features:

**The location** summarises the center of the observed values which can be interpreted as the most typical value of the variable.

**The spread** summarises how much the observed values scatter around the center.

## Location {#location}

The location of a set of observed values can be thought of as a single value that best describes the variable. It summarises where the center of the variable is located, and is sometimes called the _central tendency_. But _center_ can be interpreted in different ways, which gives rise to different measures of location. The three most common measures are the mode, the median, and the mean. ^[Here we use mean to refer to the arithmetic mean. You may also encounter the geometric mean during your studies and perhaps even the harmonic mean (although it is very rare). There are even other variations on the theme, such as weighted means.]

### The mode

The mode is simply the value that is most frequently observed. If you were to plot the values, it would refer to the peak of the plot. The problem with the mode is that it can result in several values.


```{exercise, name="The mode"}
What's the mode for nose colour in Table \@ref(tab:dataexcerpt)? What's the mode for height?
```


### The median

The median is the value that sits in the middle of the data regarding the order of magnitude. To find the median, the values of interest are arranged in ascending order and split into two equally sized groups. For example, in Table \@ref(tab:dataexcerpt), the number of antennae would be arranged as

0 | 0 | 0 | 0 | 1 | 1 | 5 | 5 | 9 | 9

We have an even number of values, producing two groups of five values each. The median is therefore assigned to be the average between the two center values, i.e. $(1+1)/2$. If we had an odd number of values, the case is easier. We just take the observation that is in the middle of the ordered sequence.

```{exercise, name="Extreme values and robustness"}
Imagine that the $10th$ value was 107 antennae, instead of 9. What would be the median now? What does the result imply for the median in terms of its robustness to outliers?
```

```{exercise, name="The Meadian"}
What's the median for nose colour in Table \@ref(tab:dataexcerpt)? What's the median for eye sight?
```


### The (arithmetic) mean

The mean is the average value of all the observations. To calculate it, we simply sum up all the values and divide it by the number of observations. The corresponding formula looks like this

$$ \bar{x} = \frac{x_1 + x_2 + \ldots + x_{n-1} + x_n}{n} = \frac{\sum_{i = 1}^n x_i}{n}, $$

where $n$ denotes the number of observations.

**Exercise**
```{exercise, name="The mean"}
Calculate the mean for number of antennae, once wih the original data and once imagining the $10th$ value was 107 antennae, instead of 9. What do the results imply for the mean in terms of its robustness to outliers? Does it make sense to calculate the mean number of antennae?
```



**Exercise**
```{exercise, name="The mean"}
What's the mean for eye sight in Table \@ref(tab:dataexcerpt)? What's the mean for nose colour?
```


These examples have shown us that the different measures of location require different mathematical operations (e.g. sorting and calculating sums). This is why not all measures of location are suitable for all variable types (cf. also Table \@ref(tab:questionstypes)).

Table: (\#tab:location-types) The types of variable determines what type of location statistic you can calculate!

| Variable Type | Mode  | Median | Mean |
|---------------|-------|--------|------| 
| Nominal       | Yes   | No     | No   |
| Ordinal       | Yes   | Yes    | No   |
| Continuous    | Yes   | Yes    | Yes  |

## Spread {#spread}

The location tells you about the typical value of a set of observations, whereas the spread tells you about how much variation there is _around_ the location. ^[Just like the location, the spread is exactly what is sounds like. Think Nutella - or schmaltz, if that's your thing - how wide can you spread it around?] Just as there are different measures of location, there are also different measures of spread.

### Proportions

Proportions are a very simple measure of spread. They are not used very often, because (unless the variable can only take two distinct values) more than a single number needs to be reported. Sometimes though, proportions are the only measure of spread that you _can_ report, see Table \@ref(tab:spread-types), below. For example, the proportions for the number of antennae are $p(0) = 0.4$ and $p(1) = p(5) = p(9) = 0.2$.

```{exercise, name="Proportions"}
What would you do if you wanted to report proportions for height?
```


### The Range

The range is defined as the distance between the minimum and the maximum of the observed values. This can either be reported as the interval or the difference between these most extreme values. As we will explain below, the minimum and maximum are sometimes referred to as $0th$ and $4th$ quartile of the data which are commonly denoted as $Q_0$ and $Q_4$, respectively. With this, the formulae for the range become $[Q_0, Q_4]$ or $Q_4 - Q_0$.

### The Interquartile Range

The interquartile range is defined as the distance between the $1st$ and the $3rd$ quartile of the data. Let's deal with quartiles first. There are five of them ^[The five quartiles $Q_0$, $Q_1$, $Q_2$, $Q_3$, $Q_4$ are sometimes also referred to as the 5-number summary.] that split the data in four even chunks. So far we have heard that the minimum is $Q_0$, and the maximum is $Q_4$. Also without noticing, we have also already learned about $Q_2$ which is the median. So there are only $Q_1$ and $Q_3$ left. For this, let's have a look at the number of antennae again. They are split into two evenly-sized groups by $Q_2$,

0 0 0 0 1 

and

1 5 5 9 9

From here it's super easy to arrive at $Q_1$ and $Q_3$ which are just the medians of the lower and upper chunk, respectively. So for the above data we have $Q_1 = 0$ and $Q_3 = 5$. Like the range, the interquartile range can come in two forms, either as the interval $[Q_1, Q_3]$ or as the difference $Q_3- Q_1$.

### The Variance

The variance is defined as an average of the squared distance of each data point from the mean^[Note the denominator which is $n-1$ instead of $n$, so the variance is \emph{almost} like taking the arithmetic mean of the squared distances $(x_i - \bar{x})^2$.]. The corresponding formula is

$$ s^2 =  \frac{\sum_{i = 1}^n (x_i - \bar{x})^2}{n-1}. $$

The values in the numerator are squared to avoid the possibility that negative and positive distances $(x_i - \bar{x})$ cancel out.
Note that the unit for the variance of a variable is the original unit squared (e.g. the variance for height would have the unit $cm^2$) which makes it difficult to work with. To simplify this, we use the standard deviation.

### The Standard Deviation

The standard deviation is the square root of the variance, i.e.
$$ s = \sqrt{s^2} $$
which provides a measurement of spread on the same scale as the original dataset. {#var-sq}

Like with the measures of location, not all measures of spread can be calculated with all variable types. Table \@ref(tab:spreadtypes) gives an overview of the possibilities.


Table: (\#tab:spreadtypes) The types of variable determines what measures of spread you can calculate.


| Variable Type | Proportions | $[Q_0, Q_4]$  or $[Q_1, Q_3]$  | $Q_4- Q_0$ or $Q_3- Q_1$ | variance or standard deviation | 
|---------------|-------------|--------------------------------|--------------------------|--------------------------------|
| Nominal       | Yes         | No                             | No                       | No                             |
| Ordinal       | Yes         | Yes                            | No                       | No                             |
| Continuous    | Yes         | Yes                            | Yes                      | Yes                            |

## Overview and outlook

We have now learned about a number of descriptive statistics and when to use them. Let's apply this to the height variable in Table \@ref(tab:questionstypes). It's displayed in Figure \@ref(fig:FakeDataForExercise).

```{r FakeDataForExercise, fig.cap = 'The values for height from 10 Martians collected at Site I.', fig.width = large_width, fig.height = 2.5, echo = FALSE}

# value <- c(188.2, 189.5, 183.0, 186.5, 191.7, 190.6, 190.9, 201.4, 191.6, 178.6)
value <- martian$Height[martian$Site == "Site I"]

# Rearrange values:
value <- value[c(1,2,3,4,5,7,6,8,10,9)]

residuals <- value - mean(value)
fake.df <- data.frame(value,
                      position = 1:10,
                      residuals)

ggplot(fake.df, aes(x = value, y = position, ymin = value - residuals, ymax = value)) +
  scale_x_continuous("Martian Height (cm)", 
                     limits = c(190,214),
                     expand = c(0,0), 
                     breaks = seq(190,214,1), 
                     labels = c(190, "", "", "", "", "",
                                196, "", "", "", "", "",
                                202, "", "", "", "", "",
                                208, "", "", "", "", "",
                                214)) +
  scale_y_continuous("Individual", limits = c(0,11.5), breaks = 1:10,  expand = c(0,0)) +
  geom_text(aes(label = as.character(value), x = value + 0.4), size = 3.0, parse = FALSE, hjust = 0) +
  geom_point(col = DarkYellow, size = 2) +
  coord_fixed(0.75) +
  theme(plot.background = element_blank(),
        panel.background = element_blank(),
        # axis.line = element_line(colour = "black"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.grid.minor.y = element_blank(),
        # panel.grid.minor.y = element_line(colour = "grey90"),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.ticks = element_line(colour = "black"),
        axis.line = element_blank())

```


```{exercise, name="Calculate values"}
Calculate $\bar{x}$, $s^2$, $s$, $Q_0$, $Q_1$, $Q_2$, $Q_3$, $Q_4$, the range and the interquartile range from Figure \@ref(fig:FakeDataForExercise)!
```

As mentioned above, both location and spread should be reported when a variable is described. The combinations that go well together are either $(\bar{x}, s)$ or $Q_2$ with the range or interquartile range.

## Outlook -- relation to statistical models

Let's talk about the combination $(\bar{x}, s)$ for a minute. If you draw $\bar{x}$ as horizontal line in Figure \@ref(fig:FakeDataForExercise) and then vertically connect each of the observed values with this line, this might remind you of a picture that you have seen before in the context of linear regression. We won't cover linear regression in this workshop until Chapter \@ref(sec:linear_regression), but let's seize the opportunity to already point out a connection between topics here!

Although they are not normally referred to as such, the vertical distances that you drew can be seen as _residuals_ that describe the discrepancy between the very simple linear model $\bar{x}$ and the data points $x_i$. The size of the individual residual $(x_i - \bar{x})$ tells you how well the model $\bar{x}$ fits the data point $x_i$, and if you were asked to provide an overall measure for the goodness of fit of the model (that avoids cancelling out between negative and positive residuals) you would quickly arrive at the standard deviation $s$!

Hold this thought, this will become an important view of the world later in the workshop!

## Boxplots

Before we move to the next chapter, we also would like to briefly talk about graphical ways of data description. A very important topic in this context is the boxplot which can be drawn using the quartile--based location/spread measures.

```{r 3-part-symm, fig.cap = 'Three visualisations of a symmetrical dataset : Height of 10 Martians from site I. _top_: A histogram, and _middle_: A dot plot, coloured according to nose colour, depicting mean and standard deviation in black. _bottom_: A box plot. The mean and median are similar, indicating that the data is symmetrically distributed around the mean.', fig.width = large_width, echo = FALSE}

# lines <- summary(martian.df$Height[martian.df$Site == "Site I"])
# lines <- lines[!names(lines) == "Mean"]
# This needs to be done manually becaus of rounding problems!
# lines <- c(176.9,   184.4,   189.0,   191.9,   203.4)

gg <- ggplot(martian[martian$Site == "Site I",], aes (x = Group, y = Height)) +
  geom_boxplot(fill = MidYellow) +
  scale_x_continuous("", limits = c(0.1,1.5), expand = c(0,0))
mData <- ggplot_build(gg)
gglines <- as.numeric(ggplot_build(gg)$data[[1]][c("ymin_final", "lower", "middle", "upper", "ymax_final")])

names(gglines) <- c("Minimum", "First Quartile", "Median", "Third Quartile", "Maximum")

#names <- rev(names)

df.x <- data.frame(gglines)

df.x <- data.frame(Statistic = row.names(df.x), Value = df.x$gglines)

# df.table <- xtable(df.x, include.rownames=FALSE, caption = "The 5-number summary depicted in a simple box-plot.", label = "FiveNumber")
# print(df.table,include.rownames=FALSE, floating.environment='margintable', latex.environments = "", table.placement = NULL)

xRangeSiteI <- c(floor(mData$data[[1]]$ymin_final),       # Get info from y axis
                 ceiling(mData$data[[1]]$ymax_final))

if (xRangeSiteI[2] %% 2) {xRangeSiteI[2] <- xRangeSiteI[2] +1 }
if (xRangeSiteI[1] %% 2) {xRangeSiteI[1] <- xRangeSiteI[1] -1 }
xRangeSiteI <- c(190, 215)
gglines <- as.numeric(ggplot_build(gg)$data[[1]][c("ymin_final", "lower", "middle", "upper", "ymax_final")])

SiteIBoxplot <- gg +
  annotate("segment",
           x = 0.1,
           xend = 1.4,
           y = gglines,
           yend = gglines,
           colour = DarkRed,
           linetype = 1) +
  annotate("text",
           x = 0.2,
           y = gglines+0.05,
           label = paste0("Q", 0:4),
           hjust=0,
           size = 2,
           angle = 0,
           colour = DarkRed) +
  scale_y_continuous("Martian Height (cm)",
                     limits = xRangeSiteI,
                     expand = c(0,0)) +
  coord_flip() +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.line = element_line(colour = "black"),
        axis.line.y = element_blank(),
        axis.text.x = element_text(colour = "black", hjust = 1),
        axis.ticks.y = element_blank(),
        axis.text.y = element_blank(),
        plot.margin = unit(c(-0.3,.2,.1,.17), "cm"),
        panel.spacing = unit(c(-0.3,0,.1,.17), "cm"))





SiteIDotplot <- ggplot(martian[martian$Site == "Site I",], aes (x = Group, y = Height, col = Nose)) +
  geom_jitter(position = position_jitter(width = .1), alpha = 0.6, size = 2.5) +
  scale_x_continuous("", limits = c(0.75,1.25), expand = c(0,0)) +
  scale_y_continuous("", limits = xRangeSiteI, expand = c(0,0)) +
  scale_colour_manual(values = c(DarkBlue, DarkRed), guide = FALSE) +
  annotate("segment",
           x = 0.75,
           xend = 1.24,
           y = gglines,
           yend = gglines,
           colour = DarkRed, linetype = 1) +
  coord_flip() +
  stat_summary(fun.y = mean, geom = "point", col = "black", size = 3.5) +
  stat_summary(fun.data = mean_sdl, geom = "errorbar", fun.args = list(mult = 1), width = 0.05, col = "black") +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        axis.line.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        plot.margin = unit(c(-.4,.2,-1,.17), "cm"), # corresponding to the top, right, bottom and left sides of the plot
        panel.spacing = unit(c(-.4,0,-1,.17), "cm"))





SiteIHistogram <- ggplot(martian[martian$Site == "Site I",], aes (x = Height, fill = Nose)) +
  geom_histogram(binwidth = 2, center = 1) +
  scale_fill_manual(values = c(MidBlue, MidRed), guide = FALSE) +
  scale_y_continuous(name = "Count", limits = c(0,3), expand = c(0,0)) +
  scale_x_continuous("Martian BMI", limits = xRangeSiteI, expand = c(0,0)) +
  annotate("segment",
           y = 0,
           yend = 3,
           x = gglines,
           xend = gglines,
           colour = DarkRed, linetype = 1) +
  theme(axis.line = element_line(colour = "black"),
        
        axis.line.x = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_line(colour = "black"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(colour = "black"),
        plot.margin = unit(c(.2,.2,-0.3,0), "cm"),
        panel.spacing = unit(c(.2,0,-0.3,0), "cm"))

ggarrange(SiteIHistogram, SiteIDotplot, SiteIBoxplot, 
          ncol = 1, nrow = 3, align = "v")

```

A box plot of the Martian height data set is shown the upper panel of Figure \@ref(fig:3-part-symm). The box encompases the IQR, and the "whiskers" extend up to, but no further than the maximum value within a _fence_ of $1.5\times$ IQR, in either direction. Data points beyond this region are represented as points. 

There are two important considerations here. First, not every graphical software program defines the default cutoff value for these points in the same way, and it can also be manually set. Second, these points _cannot automatically_ be classified outliers, although students consistently like to refer to them as such and they are even called that in many software packages. They are more accurately referred to as extreme values. Importantly, they shouldn't automatically be discarded from the analysis. Often extreme values are the interesting ones as they are unexpected and deserve a thorough investigation.

```{exercise, name="Thinking about outliers"}
If there is no universal defintion for an outlier, how would you identify and eliminate outliers?
```


Boxplots are a convenient way to visualise the distribution of a variable. For continous variables there are two other options that we'll briefly mention here, namely dot plots and histograms.

Dot plots are the most transparent visualisation method for continuous variables. Here, every data point is shown and can be coloured according to another variable. We also added the region $[\bar x - s, \bar x +s]$ to the plot. Compare the dot plot to the histogram and box plot in Figure \@ref(fig:3-part-symm).

## Skewed distributions

Ideally, we like to work with unimodal symmetrically distributed variables like Martian height shown above. Data points are centered around the mean (exhibiting only one peak) and evenly distributed on either side, becoming less frequent as we move away from the center. The mean and the median are very similar in symmetrical distributions. However, you may also observe _skewed_ distributions.

**Positively-skewed** data have a _long right tail_, with observations clustering at lower values; the mean is higher than the median. Conversly,
**Negatively-skewed** data have a _long left tail_, with observations concentrated at higher values; the mean is lower than the median.

Height at Site I is symmetrically distributed (Figure \@ref(fig:3-part-symm)). The Body-Mass Index (BMI, a combination of height and weight that provides a measure of relative size) is often positively-skewed. Below a certain value individuals are likely to die, but very high values in some individuals are possible (Figure \@ref(fig:MammoOriginalPlots)). Scores on aptitude tests, which may be too easy for the target group, are often negatively-skewed. That is, most participants are likely to have a high score (Figure \@ref(fig:neg-skew)). Another example of a negatively-skewed variable would be age at onset of a disease that appears in old-age (e.g. dementia, arthritis).

```{r MammoOriginalPlots, fig.cap = 'Three visualisation of a positively-skewed dataset : BMI of 10 Martians from Site I. _top_: A histogram, and _middle_: A dot plot, depicting mean and standard deviation, coloured according to nose colour. _bottom_: A box plot. The long right tail is evident in all visualisations and the mean is larger than the median, indicating a positive skew.', fig.width = large_width, echo = FALSE}

gg <-
  ggplot(martian[martian$Site == "Site I",], aes (x = Group, y = BMI)) +
  geom_boxplot(fill = MidYellow) +
  scale_x_continuous("", limits = c(0.1,1.5), expand = c(0,0))


mData <- ggplot_build(gg)
xRangeSiteI <-
  c(floor(mData$data[[1]]$ymin_final), # Get info from y axis
    ceiling(mData$data[[1]]$ymax_final))

if (xRangeSiteI[2] %% 2) {xRangeSiteI[2] <- xRangeSiteI[2] +1 }
if (xRangeSiteI[1] %% 2) {xRangeSiteI[1] <- xRangeSiteI[1] -1 }


gglines <- as.numeric(ggplot_build(gg)$data[[1]][c("ymin_final", "lower", "middle", "upper", "ymax_final")])

SiteIBoxplot <- gg +
  annotate("segment",
           x = 0.1,
           xend = 1.4,
           y = gglines,
           yend = gglines,
           colour = DarkRed,
           linetype = 1) +
  annotate("text",
           x = 0.2,
           y = gglines+0.05,
           label = paste0("Q", 0:4),
           hjust=0,
           size = 2,
           angle = 0,
           colour = DarkRed) +
  scale_y_continuous("Martian BMI",
                     limits = xRangeSiteI,
                     breaks = seq(xRangeSiteI[1], xRangeSiteI[2],2),
                     expand = c(0,0)) +
  coord_flip() +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.line = element_line(colour = "black"),
        axis.line.y = element_blank(),
        axis.text.x = element_text(colour = "black", hjust = 1),
        axis.ticks.y = element_blank(),
        axis.text.y = element_blank(),
        plot.margin = unit(c(-0.3,.2,.1,.17), "cm"),
        panel.spacing = unit(c(-0.3,0,.1,.17), "cm"))





SiteIDotplot <- ggplot(martian[martian$Site == "Site I",], aes (x = Group, y = BMI, col = Nose)) +
  geom_jitter(position = position_jitter(width = .1), alpha = 0.6, size = 2.5) +
  scale_x_continuous("", expand = c(0,0)) +
  scale_y_continuous("",
                     expand = c(0,0)) +
  scale_colour_manual(values = c(DarkBlue, DarkRed), guide = FALSE) +
  annotate("segment",
           x = 0.75,
           xend = 1.24,
           y = gglines,
           yend = gglines,
           colour = DarkRed, linetype = 1) +
  coord_flip(xlim = c(0.75,1.25), ylim = xRangeSiteI) +
  stat_summary(fun.y = mean, geom = "point", col = "black", size = 3.5) +
  stat_summary(fun.data = mean_sdl, geom = "errorbar", fun.args = list(mult = 1), width = 0.05, col = "black") +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        axis.line.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        plot.margin = unit(c(-.4,.2,-1,.17), "cm"),
        panel.spacing = unit(c(-.4,0,-1,.17), "cm"))





SiteIHistogram <- ggplot(martian[martian$Site == "Site I",], aes (x = BMI, fill = Nose)) +
  geom_histogram(binwidth = 2, center = 1) +
  scale_fill_manual(values = c(MidBlue, MidRed), guide = FALSE) +
  scale_y_continuous(name = "Count", limits = c(0,5), expand = c(0,0)) +
  scale_x_continuous("Martian BMI",
                     limits = xRangeSiteI,
                     expand = c(0,0)) +
  annotate("segment",
           y = 0,
           yend = 5,
           x = gglines,
           xend = gglines,
           colour = DarkRed, linetype = 1) +
  theme(axis.line = element_line(colour = "black"),
        axis.line.x = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_line(colour = "black"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(colour = "black"),
        plot.margin = unit(c(.2,.2,-0.3,0), "cm"),
        panel.spacing = unit(c(.2,0,-0.3,0), "cm"))


#Now use grid.arrange to put them all into one figure.
# grid.arrange(SiteIHistogram, SiteIDotplot, SiteIBoxplot, nrow=3)

ggarrange(SiteIHistogram, SiteIDotplot, SiteIBoxplot, 
          ncol = 1, nrow = 3, align = "v")


```


```{r neg-skew, fig.cap = 'Three visualisation of a negatively-skewed dataset: Scores on an aptitude test of 10 Martians from Site I. _top_: A histogram, and _middle_: A dot plot, depicting mean and standard deviation, coloured according to nose colour. _bottom_: A box plot. The long left tail is evident in all visualisations and the mean is lower than the median, indicating a negative skew.', fig.width = large_width, echo = FALSE}

# Aptitude 3-part

#  SiteIBoxplot
gg <-
  ggplot(martian[martian$Site == "Site I",], aes (x = Group, y = Aptitude)) +
  geom_boxplot(fill = MidYellow) +
  scale_x_continuous("", expand = c(0,0))

mData <- ggplot_build(gg)
xRangeSiteI <-
  c(floor(mData$data[[1]]$ymin_final), # Get info from y axis
    ceiling(mData$data[[1]]$ymax_final))

if (xRangeSiteI[2] %% 2) {xRangeSiteI[2] <- xRangeSiteI[2] +1 }
if (xRangeSiteI[1] %% 2) {xRangeSiteI[1] <- xRangeSiteI[1] -1 }

xRangeSiteI <- c(6,20)
gglines <- as.numeric(ggplot_build(gg)$data[[1]][c("ymin_final", "lower", "middle", "upper", "ymax_final")])

SiteIBoxplot <- gg +
  annotate("segment",
           x = 0.1,
           xend = 1.4,
           y = gglines,
           yend = gglines,
           colour = DarkRed, linetype = 1) +
  annotate("text",
           x = 0.2,
           y = gglines[1:5]+0.05, label = paste0("Q", 0:4),
           hjust=0,
           size = 2,
           angle = 0,
           colour = DarkRed) +
  scale_y_continuous("Martian score on aptitude tests",
                     limits = xRangeSiteI,
                     breaks = seq(6,20,2),
                     expand = c(0.0,0.0)) +
  coord_flip(ylim = xRangeSiteI, xlim = c(0.1,1.5)) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.line = element_line(colour = "black"),
        axis.line.y = element_blank(),
        axis.text.x = element_text(colour = "black", hjust = 1),
        axis.ticks.y = element_blank(),
        axis.text.y = element_blank(),
        plot.margin = unit(c(-0.3,.2,.1,.17), "cm"),
        panel.spacing = unit(c(-0.3,0,.1,.17), "cm"))


SiteIHistogram <- ggplot(martian[martian$Site == "Site I",], aes (x = Aptitude, fill = Nose)) +
  geom_histogram(binwidth = 2, center = 1) +
  scale_fill_manual(values = c(MidBlue, MidRed), guide = FALSE) +
  scale_y_continuous(name = "Count", expand = c(0,0)) +
  scale_x_continuous("Martian score on aptitude tests",
                     limits = xRangeSiteI,
                     expand = c(0.0,0.0),
                     breaks = xRangeSiteI) +
  coord_cartesian(xlim = xRangeSiteI, ylim = c(0,5)) +
  annotate("segment",
           y = 0,
           yend = 25,
           x = gglines,
           xend = gglines,
           colour = DarkRed,
           linetype = 1) +
  theme(
    axis.line = element_line(colour = "black"),
    axis.line.x = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_line(colour = "black"),
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
    axis.text.y = element_text(colour = "black"),
    plot.margin = unit(c(.2,.2,-0.3,0), "cm"),
    panel.spacing = unit(c(.2,0,-0.3,0), "cm"))


SiteIDotplot <-
  ggplot(martian[martian$Site == "Site I",], aes (x = Group, y = Aptitude, col = Nose)) +
  geom_jitter(position = position_jitter(width = .1), alpha = 0.6, size = 2.5) +
  # scale_x_continuous("", limits = c(0.75,1.25), expand = c(0,0)) +
  # scale_y_continuous("", limits = c(0,20), expand = c(0,0.1)) +
  scale_y_continuous("", expand = c(0,0)) +
  scale_x_continuous("", expand = c(0.0,0.0)) +
  scale_colour_manual(values = c(DarkBlue, DarkRed), guide = FALSE) +
  annotate("segment",
           x = 0.75,
           xend = 1.24,
           y = gglines,
           yend = gglines,
           colour = DarkRed,
           linetype = 1) +
  coord_flip(ylim = xRangeSiteI, xlim = c(0.75,1.25)) +
  stat_summary(fun.y = mean, geom = "point", col = "black", size = 3.5) +
  stat_summary(fun.data = mean_sdl, geom = "errorbar", fun.args = list(mult = 1), width = 0.05, col = "black" ) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        axis.line.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        plot.margin = unit(c(-0.4,.2,-1,.17), "cm"),
        panel.spacing = unit(c(-0.4,0,-1,.17), "cm"))

#Now use grid.arrange to put them all into one figure.
# library(gridExtra)
# grid.arrange(SiteIHistogram, SiteIDotplot, SiteIBoxplot, nrow=3)

ggarrange(SiteIHistogram, SiteIDotplot, SiteIBoxplot, 
          ncol = 1, nrow = 3, align = "v")

```

```{exercise, name="Skewed variables"}
Give an example of a positively or negatively-skewed variable that you will encounter in your own work.
```

Alright, now that we have an idea about how to _describe_ our _sample_, let's move onto the population.