---
title: descstat, an R Package for Computing Descriptive Statistics
date: "2025-01-08"
abstract: >
  The **descstat** package provides a set of tools to compute descriptive
  statistics, in particular, for series for which the individual
  values are recorded not as a numeric value but as belonging to 
  a numerical interval.
author:
  - name: Yves Croissant
    affiliation: LAET, Université Lumière-Lyon 2
    address:
    - MSH Lyon St-Etienne
    - 14 Avenue Berthelot
    - F-69363 Lyon Cedex 07
    - France
    url: https://laet.science
    orcid: 0000-0002-4857-7736
    email: y.croissant@univ-lyon2.fr
output: 
  html_document:
    number_sections: true
type: package
vignette: >
  %\VignetteEncoding{UTF-8}
  %\VignetteIndexEntry{desctat: an R package for descriptive statistics}
  %\VignetteEngine{knitr::rmarkdown}
---


```{r setup, echo = FALSE}
knitr::opts_chunk$set(message = FALSE, warning = FALSE,
                      fig.width = 7,
                      out.width = '80%', fig.asp = 0.6,
                      fig.align = "center")
options(
    htmltools.dir.version = FALSE, formatR.indent = 2, width = 55, digits = 4,
    tibble.print_max = 5, tibble.print_min = 5)
```

**R** offers many tools to analyze the univariate or bivariate
distribution of series. This includes `table` and `prop.table` in base
**R**, `group_by/summarise` and `count` in **dplyr**. However,
these functions are somehow frustrating as some very common tasks,
like:

- adding a total,
- computing relative frequencies or densities instead of counts,
- merging high values of integer series in a `>=N` category,
- computing cumulative distributions, 
- computing descriptive statistics,

are tedious. Moreover, to our knowledge, **R** offers weak support for
numerical series for which the numerical value is not known at the
individual level, but only the fact that this value belongs to an
interval. **descstat** is a small dependency-free package that is
intended to provide user-friendly tools to perform this kind of
operations. More specifically, **descstat** provides\:

- a `bin` class,
- a `freq_table` function to construct a frequency table, and some
  methods to plot and to compute descriptive statistics,
- a `cont_table` function to compute a contingency table for two
  series.

These function are writen in the **tidyverse** style, which means that
 series can be selected without quotes (but note that characters can
 also be used). Plots are computed using the **tinyplot** package and
 the **descstat** package provides several **tinyplot**'s types. We
 therefore load both packages and we use **tinyplot**'s `"clean"`
 theme thorough.

```{r load_pkgs}
library(descstat)
library(tinyplot)
tinyplot::tinytheme("clean")
```

# Bins

We'll call bin a series that has a limited set of distinct values. A
bin is created using the `bin` function. The `bin` function returns a
series of class `bin`. The input can be either :

- a continuous series,
- a categorical series,
- an integer series,
- a continuous bin series, i.e. a series where a continuous numerical
  variable is coded as numerical intervals.


## Continuous numerical series

When a continuous numerical series is provided, the `bin` function
coerces it to a continuous bin series, using the `base::cut` function,
resulting in a factor, with labels having a strict format, like
`[10,20)`, `(20,30]`, `[50,Inf)`.

<!-- - an oppening bracket that should be either `(` or `[` for -->
<!--   respectively a bin open or closed on the left, -->
<!-- - a first numerical value (the lower bound of the bin), -->
<!-- - a comma, -->
<!-- - a second numerical value (the upper bound of the bin), -->
<!-- - a closing bracket that should be either `)` or `]` for -->
<!--   respectively a bin open or closed on the right. -->
  
As `base::cut`, the `bin` function has a `breaks` argument (a numeric
containing the breaks) and a `right` argument (if `TRUE`, the default,
bins are closed on the right).

<!-- ```{r using_cut} -->
<!-- z <- c(1, 5, 10, 12, 4, 9, 8) -->
<!-- bin1 <- cut(z, breaks = c(1, 8, 12), right = FALSE) -->
<!-- bin2 <- cut(z, breaks = c(1, 8, 12), right = TRUE) -->
<!-- bin3 <- cut(z, breaks = c(1, 8, 12, Inf), right = FALSE) -->
<!-- data.frame(z, bin1, bin2, bin3) -->
<!-- ``` -->

<!-- Note that\ : -->

<!-- - the value of 8 is included in the `[8,12)` bin when `right = FALSE` -->
<!--   and in the `(1,8]` bin when `right = TRUE`, -->
<!-- - the value of 1 results in a `NA` when `right = TRUE` as the lowest -->
<!--   value of `breaks` is 1, -->
<!-- - the value of 12 results in a `NA` when `right = FALSE` as 12 is the -->
<!--   highest value of `breaks`. -->

The `breaks` argument of the `bin` function is not mandatory. If it is
missing, a vector of breaks is auto-processed. Using the `padova` data
set, which contains the price and the characteristics of a sample of
housings in Padova, we get:

```{r }
#| collapse: true
padova$price |> range()
bin(padova$price) |> head()
```

If `breaks` is provided, its first (last) values can be
either:

- outside the range of the series; in this case the lower bound of the
  first bin and the upper bound of the last bin are these values,
- inside the range of the series; in this case the lower bound of the
  first bin is `0` (and not the lower value of the series) and the
  upper bound of the last bin is `Inf`,
  
as illustrated in the two following examples:  

```{r freq_table_continuous}
bin(padova$price, breaks = c(250, 500, 750)) |> head()
bin(padova$price, breaks = c(30, 250, 500, 750, 1000)) |> head()
```

## Continuous bin series

When the argument of `bin` is already a character or a factor
containing numerical intervals, the correctness of the values are
checked and replaced by `NA`s in case of problems and a factor is
returned with the levels in the relevant order.  The initial number of
intervals can be reduced using the `breaks` argument. Consider for
example the `wage` series of the `wages` data set:

```{r }
head(wages$wage)
```
Providing a new vector of breaks, we get:

```{r }
wages$wage |> bin(breaks = c(1, 4, 8)) |> head()
```

If a unique break is provided all the values greater or equal to
`breaks` are collapsed in a unique interval:

```{r }
wages$wage |> bin(breaks = 8) |> head()
```

## Integer series

A series is considered to be an integer series if two conditions are
met:

- the first one is that it should "looks" like an integer, i.e. the
  difference between `x` and `round(x)` should be lower than a given
  level,
- the number of values should be reasonably low: for example, the
  default behaviour is to consider that a series is an integer series
  if the number of distinct values is lower than 50.
  
The `bin` function doesn't do anything in this case except if the
`max` argument is set: in this case, all the values greater or equal
than this value are coerced to a unique value which is the mean of the
variable in this interval. Using the `children` series of the `rgp`
data set and setting the `max` argument to 3, we get:

```{r }
rgp$children |> head(12)
rgp$children |> bin(max = 3) |> head(12)
```

## Categorical series

If the input is a series which is either a character or a factor and
doesn't follow the strict syntax of continuous bins, it is considered
to be a categorical variable and is coerced if necessary to a
factor. Consider for example the `sector` series of the `wages` data
set:

```{r }
wages$sector |> head(5)
```

Without any further argument, the `bin` function just returns the
series:

```{r results = "hide"}
wages$sector |> bin()
```

The `levels` argument can be used to reorder, rename or collapse some
of the initial levels. It is a (potentially named) character,
except for collapsing where it should be a named list:

To get a different order of the levels:

```{r }
wages$sector |> bin(levels = c("administration", "business", "services", 
                             "industry", "building")) |> head()

```

To rename some of the levels:

```{r }
wages$sector |> bin(levels = c(government = "administration", "business", "services", 
                             "industry", construction = "building")) |> head()

```
To collapse levels:

```{r }
wages$sector |> 
    bin(levels = list(pub = "administration", priva = c("industry", "building"),
                      privb = c("business", "services"))) |>  head()
```

## From continuous bins to numerical values

The real strength of the `bin` class for continuous bin series is to
allow calculus on the underlying numerical values. To perform this
task, an `as_numeric` function is provided which returns a numeric
series. The use of the `as_numeric` function is illustrated using the
first 6 values of `size` series of the `wages`:


```{r wsize}
wsize <- head(wages$size)
wsize
```

Using `as_numeric` without further argument simply returns the lower
bound of the intervals\:

```{r as_numeric_base}
wsize |> as_numeric()
```

but a `pos` argument can be used to return any value in the range of
the bin, by defining the relative position, i.e.:

- `pos = 0` (the default) for the lower bound,
- `pos = 1` for the upper bound,
- `pos = 0.5` for the center of the bin.

```{r as_numeric_ex}
wsize |> as_numeric(pos = 1)
wsize |> as_numeric(pos = 0.5)
```

While computing some statistics on bins, it is customary to consider
that, for all the individuals belonging to a specific interval, the
value is just the center (i.e., all the individual values are equal to
15 for individuals belonging to the `[10,20)` interval), or,
equivalently, in terms of mean, that the values are uniformly
distributed in the interval. This value is returned by `as_numeric`
with `pos = 0.5`, but there is a specific problem for the last
interval if, as it is the case here, it is open to infinity on the
right. In this case, the default behavior of `as_numeric` is to set
the width of the last interval to the width of the penultimate one. As
the width of the `[100,250)` interval is 150, the width of the last one is
also set to 150, which mean an upper bound of 400 and a center value
of 325. This behavior can be changed either by setting the `wlast`
argument to a value different than 1, which is interpreted as a
multiple of the width of the penultimate interval. For example, `wlast
= 2` means that the width of the last interval is set to 2 times the
one of the penultimate interval, which means 300, and the resulting
interval is `[250,550)` with a center value of 400.

```{r as_numeric_wlast}
wsize |> as_numeric(pos = 0.5, wlast = 2)
```
The same result can be obtained by directly setting the center of the
last bin using the `xlast` argument\:

```{r as_numeric_xlast, results = 'hide'}
wsize |> as_numeric(pos = 0.5, xlast = 400)
```

Finally, a specific numerical value for the first bin can also be set using
the `xfirst` argument\:

```{r as_numeric_xfirst}
wsize |> as_numeric(pos = 0.5, xlast = 400, xfirst = 2)
```

`xfirst`, `xlast` and `wlast` are also arguments that can be used in
`bin`. In this case, they are stored as attributes of the resulting
`bin` series and are used while calling `as_numeric`.

```{r bin_attributes}
wsize |> bin(xlast = 400, xfirst = 2) |> as_numeric(pos = 0.5)
```

# Computing frequency tables

## Discrete series

A frequency table summarizes the univariate distribution of a `bin`
series. This task can be performed using `base:table` or, with
**tidyverse**, using `dplyr::count`. For example, the `rgp` data set,
which is an extract of the French census, contains the number of
children in households:

```{r rgp}
rgp |> head(5)
```

Using `table`, coercing to a data frame and setting relevant names, we
get:

```{r rgp_count}
table(rgp$children) |> as.data.frame() |> setNames(c("children", "n"))
```

The `descstat::freq_table` function performs the same task and returns
by default exactly the same data frame as previously, except that the
resulting `children` series is of class `bin`:

```{r freq_table_base, results = 'hide'}
rgp |> freq_table(children)
```

Several further arguments are provided which can improve the result:

- `f` is a character containing one or several letters and indicates
  what kind of frequencies should be computed\:
  
    - `n` for counts or absolute frequencies (the default),
    - `f` for the (relative) frequencies,
    - `p` for the percentage (i.e., $f\times 100$),
	- `N`, `F` and `P` for the cumulative values of `n`, `f` and
      `p`.
- `total` a boolean, if `TRUE`, a total is returned,
- `max`, suitable for an integer series only, is an integer which is the
  last value presented in the frequency table, e.g., `max = 3` creates
  a last line which is `>=3`.
  
The following command uses the `f` argument with all the possible letters.

```{r freq_table_freqs}
rgp |> freq_table(children, "nfpNFP")
```

As there are few occurrences of families with more than 3 children, we
set `max = 3` and add a total by setting `total` to `TRUE`.

```{r freq_table_max}
rgp |> freq_table(children, max = 3, total = TRUE)
```

Note that the `children` series is still numeric: all the values
greater or equal to 3 are replaced by the mean value in this interval
($3.344$) and the value of the series for the total is `NA`. To get a
printable version of the table, the `print` argument of `freq_table`
should be set to `TRUE`:

```{r freq_table_max_print}
rgp |> freq_table(children, max = 3, total = TRUE, print = TRUE)
```

## Continuous series

Frequency tables can also be computed for continuous series, provided
either as intervals or as numerical values.  Applying
`descstat::freq_table` to the `size` series of the `wages` data set,
we get\:

```{r freq_table_size}
wages |> freq_table(size)
```

A `breaks` argument can be provided which is then passed internally to
the `bin` function (see section 1.2) in order to create a bin series
with less modalities :

```{r freq_table_breaks}
wages |> freq_table(size, breaks = c(20, 250))
```

which is equivalent to

```{r freq_table_breaks2, results = 'hide'}
wages |>
    transform(size = bin(size, breaks = c(20, 250))) |>
    freq_table(size)
```

A frequency table with numerical bins can also be created from a
continuous numerical series, as `price` in the `padova` data set \:
<!-- (see section 1.1 for the use of the `breaks` argument) \: -->

```{r range_padova, collapse = TRUE}
padova |> freq_table(price)
```

The `f` argument may contain further letters than for the discrete
series case\:

- `d` for density, which is the relative frequency divided by the
  width of the interval,
- `m` for the mass frequency (the mass is the relative frequency of the
 mass of the variable, i.e., of the product of the frequency and the
 value of the variable),
- `M` for cumulative mass frequencies.

```{r freq_table_bins}
wages |> freq_table(size, "dmM", breaks = c(20, 100, 250))
```

Other values of the variable can be included in the table using the
`vals` argument, which is a character including some of the following
letters:

- `l` for the lower bound of the interval,
- `u` for the upper bound of the interval,
- `a` for the width of the interval.

```{r freq_table_bins_vals}
wages |> freq_table(size, "p", vals = "xlua", breaks = c(20, 100, 250), wlast = 2)
```

## Dealing with frequency tables

`freq_table` is not only designed for data sets containing individual
observations but also for frequency tables, like the `income` data
set\:

```{r income}
income |> head(5)
```

which presents the distribution of income in France with 25 intervals;
`inc_class` contains ranges of yearly income (in thousands of euros),
`number` is the number of households (in millions) and `tot_inc` is
the mass of income in each bin (in thousands of million). To use
`freq_table` to read this frequency table, a further argument `freq`
is required and should indicate which column of the data frame contains
the frequencies\:

```{r freq_table_freq}
income |> freq_table(inc_class, freq = number,
                     breaks = c(50, 100, 500, 1000))
```
If one column of the data frame contains the mass of the variables (the
`tot_inc` series in the `income` data set), it can be indicated using
the `mass` argument\:

```{r freq_table_mass}
income |> freq_table(inc_class, freq = number,
                     mass = tot_inc,
                     breaks = c(50, 100, 500, 1000))
```

The center of the bins are then calculated by dividing the mass by
the frequency, which result, by definition, to the exact mean value for
every bin.

## Dealing with weights

The `weights` argument can be used to mimic the population when the
data set contains a vector of weights for every observations. For
example, the `employment` table contains a series of weights called
`weights`. Using `freq_table` for the `age` series and setting the
`weights` argument to the series of weights, we get:

```{r cont_table_weights}
employment |> freq_table(age, weights = weights, total = TRUE, print = TRUE)
```

and  the resulting `n` series is the sum of the weights for
every modality of the series.


## Univariate descriptive statistics

Descriptive statistics can easily be computed applying functions to
`freq_table` objects. The problem is that only a few statistical
functions of **R**'s **base** and **stats** packages are generic
(`mean`, `median` and `quantile`). For these, methods where written
for `freq_table` objects. For the other ones, we had to create new
functions. Table 1 indicates the **R** functions
and the corresponding **descstat** functions.

```{r functions, echo = FALSE}
data.frame(R        = c("mean", "median", "quantile", "var",
                        "sd",    "mad", "", "", "", "", ""),
           descstat = c("mean", "median", "quantile", "variance",
                        "stdev", "madev", "modval", "medial", "gini",
                        "skewness", "kurtosis")) |> 
    knitr::kable(caption = "Table 1: Functions for descriptive statistics",
                 booktabs = TRUE)
```

To compute the central values statistics of the distribution of
`wages` we use:

```{r central_stat, collapse = TRUE}
z <- wages |> freq_table(wage)
z |> mean()
z |> median()
z |> modval()
```

`median` returns a value computed using a linear interpolation and
`modval` returns the mode, which is a one line data frame containing
the interval, the center of the interval and the frequency for the
mode.

For the dispersion statistics,^[`madev` computes the mean absolute
variation.], we get:

```{r disp_stat, collapse = TRUE}
z |> stdev()
z |> variance()
z |> madev()
```
For the quantiles, the argument `y` can be used to compute the
quantiles using the values of the variable (`y = "value"`, the
default) or the masses (`y = "mass"`):

```{r quantiles, collapse = TRUE}
z |> quantile(probs = c(0.25, 0.5, 0.75))
z |> quantile(y = "mass", probs = c(0.25, 0.5, 0.75))
```

The quantile of level 0.5 is the median in the first case, the medial
in the second case\:

```{r median_medial, collapse = TRUE}
z |> median()
z |> medial()
```

`gini` computes the Gini coefficient of the series:

<!-- [Twice the grey area in figure \@ref(fig:lorenz).] -->

```{r gini_coef, collapse = TRUE}
z |> gini()
```

`skewness` and `kurtosis` compute Fisher's shape statistics\:

```{r shape_stat, collapse = TRUE}
z |> skewness()
z |> kurtosis()
```

All these functions also work for counts. Of course, for categorical
series, all these functions are irrelevant except the one which
computes the mode\:

```{r modval_wages}
wages |> freq_table(sector) |> modval()
```

# Ploting a frequency table

## Discrete series

For an integer or a categorical series, the most natural way to plot
the distribution is to use a bar plot. For categorical series, this
can be done using the built-in `barplot` type of **tinyplot**, either
using a pre-computed frequency table (see figure 1),

```{r barplot, fig.cap = "Figure 1:  Bar plot for a frequency table"}
s <- freq_table(wages, sector, f = "f")
tinyplot(f ~ sector, data = s, type = "barplot", flip = TRUE, yaxl = "%", ylab = "Frequencies")
```

or using a raw table:
<!-- (see figure \@ref(fig:barplotraw)). -->

```{r barplotraw, fig.cap = "Bar plot for a raw table", eval = FALSE}
tinyplot(~ sector, data = wages, type = "barplot", flip = TRUE, ylab = "Frequencies")
```

For integer series, if the `max` argument has been used, the `print`
argument should be set to `TRUE`, so that the value is of the form `">=3"` and
not `3.34` (see figure 2).

```{r barplotint, fig.cap = "Figure 2: Bar plot for integer series"}
i <- freq_table(rgp, children, max = 3, f = "f", print = TRUE)
tinyplot(f ~ children, data = i, type = "barplot", flip = TRUE, yaxl = "%")
```

Tinyplot's types provided by **descstat** follow the same logic. They
have a formula/data interface that can be either:

- a one-side formula and a raw data set,
- a two-side formula and a frequency table.

An (unadvised) alternative to bar plots are pie charts. They can be
produced using `descstat::type_pie()` (see figure 3).


```{r pie2, fig.asp = 1, out.width = "50%", fig.cap = "Figure 3: A pie chart"}
tinyplot(~ sector, wages, type = type_pie())
```
or:

```{r pie2bis, fig.asp = 1, out.width = "50%", eval = FALSE}
ws <- freq_table(wages, sector)
tinyplot(n ~ sector, ws, type = type_pie())
```


The position of the labels can be set using the `position` argument
(the default value is 0.5), the initial angle using the `init.angle`
and the size of the foreground using the `radius` argument which is
useful when there is insufficient place to write some labels (see
figure 4).^[These arguments are the same as
`graphics::pie`.]

```{r pie3, fig.asp = 1, out.width = "50%", fig.cap = "Figure 4: A pie chart with customized initial angle and labels positions"}
tinyplot(~ sector, wages,
         type = type_pie(position = c(1.1, 0.7), init.angle = 10, f = "p",
                         radius = 1.3, pal = "Tableau 10"))
```
A donut chart is obtained by setting the `hole` argument to a value
strictly between 0 and 1 (see figure 5).

```{r donut, fig.asp = 1, out.width = "50%", fig.cap = "Figure 5: Donut chart"}
tinyplot(~ sector, wages,
         type = type_pie(position = 1.1, init.angle = 10,
                         radius = 1.3, hole = 0.3, pal = "Set3"))
```

The cumulative distribution can be plotted using the
`descstat::type_cumul` type (see figure 6). The
`expand` argument is a fraction of the range of the series and set the
width of the segments for which the cumulative distribution is 0 or 1.

```{r cummulative, fig.cap = "Figure 6: Cummulative distribution"}
tinyplot(~ children, data = rgp, type = type_cumul(expand = 0.1, max = 3))
```

Using a frequency table, we get:

```{r cummulativebis, eval = FALSE}
rc <- freq_table(rgp, children, f = "F", max = 3)
tinyplot(F ~ children, data = rc, type = type_cumul(expand = 0.1))
```

## Continuous series

Relevant plots for continuous series are very different from those
suitable for discrete or categorical series. The two most popular
plots are histogram and frequency polygon, respectively obtained by
setting the `type` argument of `tinyplot` to `descstat::type_histo`
(see figure 7),

```{r histogram, fig.cap = "Figure 7: Histogram"}
tinyplot(~ wage, data = wages,
         type = type_histo(breaks = c(10, 20, 30, 40, 50)))
```

and `descstat::type_freqpoly` (see figure 8).

```{r freqpoly, fig.cap = "Figure 8: Frequency polygons"}
tinyplot(~ wage, data = wages,
         type = type_freqpoly(breaks = c(10, 20, 30, 40, 50)))
```
Using a frequency table, the same plots can be obtained using:


```{r freqpolybis, eval = FALSE}
ww <- freq_table(wages, wage, f = "d",
                 breaks = c(10, 20, 30, 40, 50))
tinyplot(d ~ wage, data = ww, type = type_freqpoly())
```

As for integer series (see figure 6), the
cummulative distribution can be plotted using the
`descstat::` `type_cumul` type (see figure 9).

```{r cumm2, fig.cap = "Figure 9: Cummulative distribution"}
tinyplot(~ wage, data = wages, type = type_cumul(expand = 0.05))
```

Another popular plot for continuous series is the Lorenz curve, which
indicates the relation between the cumulative distributions of the
frequencies and of the masses of the series . A Lorenz curve is
obtained by setting the `type` argument to `descstat::type_lorenz`
(see figure 10).

```{r lorenz, fig.cap = "Figure 10: Lorenz curve"}
tinyplot(~ wage, wages, type = type_lorenz(), col = "lightgrey")
```

or, using a frequency table:

```{r lorenzbis, eval = FALSE}
ww <- freq_table(wages, wage, f = "FM")
tinyplot(M~ F, ww, type = type_lorenz(), col = "lightgrey")
```

# Contingency tables

With base **R**, a contingency table can be computed using the `table`
function with two categorical series.  To illustrate the
computation of contingency tables, we'll use the `size` and `wage`
series of the `wages` table and we'll first reduce the number of
intervals:

```{r wages2}
wages2 <- wages |>
    transform(size = bin(size, breaks = c(20, 50, 100)),
              wage = bin(wage, breaks = c(10, 30, 50)))
```

Using the `table` function, we get a contingency table in "wide"
format:

```{r }
with(wages2, table(wage, size))
```

The result is not a data frame, but a matrix. Coercing it to a data
frame, we get a contingency table in "long" format:

```{r }
with(wages2, table(wage, size)) |> as.data.frame() |> head()
```

## Creating contingency tables using `cont_table`

The same contingency table can be obtained using
`descstat::cont_table`:

```{r cont_table}
wages2 |> cont_table(wage, size) |> head()
```

The `descstat::wide` function transforms the contingency table in wide
format:

```{r cont_table_wide}
wages2 |> cont_table(wage, size) |> wide()
```

as for `freq_table`, central values for the first and the last class
can be set using arguments `xfirst#`, `xlast#` and `wlast#`, where `#`
is equal to 1 or 2 for the first and the second series indicated in
the `cont_table` function.

## Plotting a contingency table

A contingency table can be plotted using the
`descstat::type_cont.table` type of **tinyplot** (see figure
11). This is simply a scatterplot, with the size
of the points related to the frequencies. The `size` argument
indicates the size of the smallest and of the largest point.


```{r conttableplot, fig.cap = "Figure 11: Contingency table plot"}
tinyplot(wage ~ size, data = wages2,
         type = type_cont.table(size = c(0.5, 5)),
         pch = 16)
```

## Computing the distributions from a contingency table

Denote $n_{ij}$ the count of the cell corresponding to the $i$^th^
modality of the first variable and the $j$^th^ modality of the second
one. Three distributions can be computed:

- the joint distribution is obtained by dividing $n_{ij}$ by the
  sample size,
- the marginal distribution is obtained by summing the counts row-wise
  $n_{i.}=\sum_{j}n_{ij}$ for the first variable and column-wise
  $n_{.j} = \sum_{i}n_{ij}$ for the second one, and dividing them by
  the sample size,
- the conditional distribution is obtained by dividing the joint by
  the marginal distribution.
  
The `joint`, `marginal` and `conditional` functions return these three
distributions. The last two require an argument `y` which is one of the
two series of the `cont_table` object.

```{r distributions}
wht <- wages2 |> cont_table(size, wage)
wht |> joint() |> wide()
wht |> marginal(size)
wht |> conditional(size) |> wide()
```

Note that `marginal`, as it returns an univariate distribution, is
coerced to a `freq_table` object.

## Computing descriptive statistics

Descriptive statistics can be computed using any of the three
distributions. Using the joint distribution, we get a data frame
containing two columns for the two series.

```{r distributions_stat, collapse = TRUE}
wht |> joint() |> mean()
wht |> joint() |> stdev()
wht |> joint() |> variance()
wht |> joint() |> modval()
```
The same (univariate) statistics can be obtained using the marginal
distribution:

```{r marginal_mean, collapse = TRUE}
wht |> marginal(size) |> mean()
```
or even more simply considering the univariate distribution computed
by `freq_table`:

```{r marginal_mean_freq_table, collapse = TRUE}
wages2 |> freq_table(size) |> mean()
```

The `mean`, `stdev` and `variance` methods are actually only usefull
when applied to a conditional distribution; in this case, considering
for example the conditional distribution of the `wage` variable, there
are as many values returned than the number of modalities of the
second (conditioning) variable:

```{r conditional_stat}
wht |> conditional(wage) |> mean()
wht |> conditional(wage) |> variance()
```

## Analysis of variance


The total variance of $X$ can be written as the sum of:

- the residual variance, *ie* the mean of the conditional variances.
- the explained variance, *ie* the variance of the conditional means,

Denoting $\bar{x}_j$ and $s^2_{x_j}$ the conditional mean and the
conditional variance and $f_{.j}$ the marginal frequencies:

$$
s_{x}^2 = \sum_j f_{.j} s^2_{x_j} + \sum_j f_{.j} (\bar{x}_j -
\bar{\bar{x}}) ^ 2
$$

The decomposition of the variance can be computed by joining tables
containing the conditional moments and the marginal distribution of
the conditioning variable and then applying the formula:

```{r anova_computations}
cm <- wht |> conditional(wage) |> mean()
cv <- wht |> conditional(wage) |> variance()
md <- wht |> marginal(size)
z <- merge(md, cm) |> merge(cv)
om <- sum(z$f * z$mean)
ev <- sum(z$f * (z$mean - om) ^ 2)
rv <- sum(z$f * z$variance)
tv <- ev + rv
c(om = om, ev = ev, rv = rv, tv = tv)
```

Or more simply using the `anova` method for `cont_table` objects:

```{r anova}
wht_wage <- wht |> anova(wage)
wht_wage
```

which has a `summary` method which computes the different elements of
the decomposition:

```{r anova_summary}
wht_wage |> summary()
```
and especially the correlation ratio, which is obtained by dividing
the explained variance by the total variance. 

The anova plot of `wage` on `size` is obtained using `descstat::type_anova`
(see figure 12).

```{r anovaplot, fig.cap = "Figure 12: Anova plot"}
tinyplot(wage ~ size, data = wages2, type = type_anova())
```
The points are the mean values of `wage` for every value of `size` and
the error bars represents the confidence interval for a given level,
that can be set using the `level` argument (the default value is 0.95).


## Linear regression


For the joint distribution, two other functions are provided,
`covariance` and `correlation` (which are the equivalent of the
non-generic `stats::cov` and `stats::cor` functions) to compute the
covariance and the linear coefficient of correlation.

```{r contingent_covariance, collapse = TRUE}
wht |> joint() |> covariance()
wht |> joint() |> correlation()
```

The regression line can be computed using the `regline` function:

```{r regline, collapse = TRUE}
rl <- regline(wage ~ size, wht)
rl
```
which returns the intercept and the slope of the regression of `wage`
on `size`. We can then add the regression line to figure
11 (see figure 13).


```{r reglineplot, fig.cap = "Figure 13: Regression line"}
tinyplot(wage ~ size, data = wages2, type = type_cont.table(size = c(0.5, 5)),
         legend = FALSE, pch = 16)
a <- cont_table(wages2, wage, size)
cfs <- regline(wage ~ size, data = a)
tinyplot_add(type = type_abline(a = cfs[1], b = cfs[2]), col = "blue")
```
