--- 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)`. 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). 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 \: ```{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: ```{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: ```{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") ```