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:

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:

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.

library(descstat)
library(tinyplot)
tinyplot::tinytheme("clean")

1 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 :

1.1 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:

padova$price |> range()
## [1]  35 950
bin(padova$price) |> head()
## [1] (0,100]   (200,300] (200,300] (100,200] (100,200]
## [6] (400,500]
## 6 Levels: (0,100] (100,200] (200,300] ... (500,Inf]

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:

bin(padova$price, breaks = c(250, 500, 750)) |> head()
## [1] (0,250]   (0,250]   (0,250]   (0,250]   (0,250]  
## [6] (250,500]
## Levels: (0,250] (250,500] (500,750] (750,Inf]
bin(padova$price, breaks = c(30, 250, 500, 750, 1000)) |> head()
## [1] (30,250]  (30,250]  (30,250]  (30,250]  (30,250] 
## [6] (250,500]
## 4 Levels: (30,250] (250,500] ... (750,1e+03]

1.2 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 NAs 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:

head(wages$wage)
## [1] [14,16)  [4,6)    [40,50)  [8,10)   [50,Inf)
## [6] [14,16) 
## 23 Levels: [0.2,0.5) [0.5,1) [1,1.5) [1.5,2) ... [50,Inf)

Providing a new vector of breaks, we get:

wages$wage |> bin(breaks = c(1, 4, 8)) |> head()
## [1] [8,Inf) [4,8)   [8,Inf) [8,Inf) [8,Inf) [8,Inf)
## Levels: [0,1) [1,4) [4,8) [8,Inf)

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

wages$wage |> bin(breaks = 8) |> head()
## [1] [8,Inf) [4,6)   [8,Inf) [8,Inf) [8,Inf) [8,Inf)
## 9 Levels: [0.2,0.5) [0.5,1) [1,1.5) [1.5,2) ... [8,Inf)

1.3 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:

rgp$children |> head(12)
##  [1] 0 0 0 1 0 0 2 1 2 1 1 5
rgp$children |> bin(max = 3) |> head(12)
##  [1] 0.000 0.000 0.000 1.000 0.000 0.000 2.000 1.000
##  [9] 2.000 1.000 1.000 3.344

1.4 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:

wages$sector |> head(5)
## [1] industry       administration business      
## [4] administration services      
## 5 Levels: industry building business ... administration

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

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:

wages$sector |> bin(levels = c("administration", "business", "services", 
                             "industry", "building")) |> head()
## [1] industry       administration business      
## [4] administration services       business      
## 5 Levels: administration business ... building

To rename some of the levels:

wages$sector |> bin(levels = c(government = "administration", "business", "services", 
                             "industry", construction = "building")) |> head()
## [1] industry   government business   government
## [5] services   business  
## 5 Levels: government business services ... construction

To collapse levels:

wages$sector |> 
    bin(levels = list(pub = "administration", priva = c("industry", "building"),
                      privb = c("business", "services"))) |>  head()
## [1] priva pub   privb pub   privb privb
## Levels: pub priva privb

1.5 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:

wsize <- head(wages$size)
wsize
## [1] [1,10)    [50,100)  [20,50)   [250,Inf) [250,Inf)
## [6] [20,50)  
## 6 Levels: [1,10) [10,20) [20,50) ... [250,Inf)

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

wsize |> as_numeric()
## [1]   1  50  20 250 250  20

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.
wsize |> as_numeric(pos = 1)
## [1]  10 100  50 400 400  50
wsize |> as_numeric(pos = 0.5)
## [1]   5.5  75.0  35.0 325.0 325.0  35.0

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.

wsize |> as_numeric(pos = 0.5, wlast = 2)
## [1]   5.5  75.0  35.0 400.0 400.0  35.0

The same result can be obtained by directly setting the center of the last bin using the xlast argument:

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:

wsize |> as_numeric(pos = 0.5, xlast = 400, xfirst = 2)
## [1]   2  75  35 400 400  35

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.

wsize |> bin(xlast = 400, xfirst = 2) |> as_numeric(pos = 0.5)
## [1]   2  75  35 400 400  35

2 Computing frequency tables

2.1 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:

rgp |> head(5)
##   cars rooms children   type
## 1    2     4        0 couple
## 2    2     8        0 couple
## 3    1     3        0 couple
## 4    1     3        1 couple
## 5    2     2        0 couple

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

table(rgp$children) |> as.data.frame() |> setNames(c("children", "n"))
##   children   n
## 1        0 317
## 2        1 272
## 3        2 260
## 4        3 112
## 5        4  31
## 6        5   6
## 7        6   1
## 8        9   1

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:

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.

rgp |> freq_table(children, "nfpNFP")
##   children   n     f    p    N     F     P
## 1        0 317 0.317 31.7  317 0.317  31.7
## 2        1 272 0.272 27.2  589 0.589  58.9
## 3        2 260 0.260 26.0  849 0.849  84.9
## 4        3 112 0.112 11.2  961 0.961  96.1
## 5        4  31 0.031  3.1  992 0.992  99.2
## 6        5   6 0.006  0.6  998 0.998  99.8
## 7        6   1 0.001  0.1  999 0.999  99.9
## 8        9   1 0.001  0.1 1000 1.000 100.0

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.

rgp |> freq_table(children, max = 3, total = TRUE)
##   children    n
## 1    0.000  317
## 2    1.000  272
## 3    2.000  260
## 4    3.344  151
## 5       NA 1000

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:

rgp |> freq_table(children, max = 3, total = TRUE, print = TRUE)
##   children    n
## 1        0  317
## 2        1  272
## 3        2  260
## 4      >=3  151
## 5    Total 1000

2.2 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:

wages |> freq_table(size)
##        size     x   n
## 1    [1,10)   5.5 207
## 2   [10,20)  15.0  90
## 3   [20,50)  35.0 124
## 4  [50,100)  75.0 113
## 5 [100,250) 175.0 124
## 6 [250,Inf) 325.0 342

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 :

wages |> freq_table(size, breaks = c(20, 250))
##        size   x   n
## 1    [0,20)  10 297
## 2  [20,250) 135 361
## 3 [250,Inf) 365 342

which is equivalent to

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 : 

padova |> freq_table(price)
##       price   x   n
## 1   (0,100]  50  47
## 2 (100,200] 150 469
## 3 (200,300] 250 230
## 4 (300,400] 350 150
## 5 (400,500] 450  87
## 6 (500,Inf] 550  59

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.
wages |> freq_table(size, "dmM", breaks = c(20, 100, 250))
##        size   x         d       m       M
## 1    [0,20)  10 0.0148500 0.01979 0.01979
## 2  [20,100)  60 0.0029625 0.09477 0.11457
## 3 [100,250) 175 0.0008267 0.14463 0.25920
## 4 [250,Inf) 325 0.0022800 0.74080 1.00000

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.
wages |> freq_table(size, "p", vals = "xlua", breaks = c(20, 100, 250), wlast = 2)
##        size   x    p   l   u   a
## 1    [0,20)  10 29.7   0  20  20
## 2  [20,100)  60 23.7  20 100  80
## 3 [100,250) 175 12.4 100 250 150
## 4 [250,Inf) 400 34.2 250 550 300

2.3 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:

income |> head(5)
##   inc_class number tot_inc
## 1    [0,10)  8.763   36.33
## 2   [10,12)  2.108   23.25
## 3   [12,15)  3.396   46.24
## 4   [15,20)  6.000  104.44
## 5   [20,30)  6.979  171.52

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:

income |> freq_table(inc_class, freq = number,
                     breaks = c(50, 100, 500, 1000))
##     inc_class    x         n
## 1      [0,50)   25 34.015929
## 2    [50,100)   75  3.505522
## 3   [100,500)  300  0.785294
## 4 [500,1e+03)  750  0.018466
## 5 [1e+03,Inf) 1250  0.007766

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:

income |> freq_table(inc_class, freq = number,
                     mass = tot_inc,
                     breaks = c(50, 100, 500, 1000))
##     inc_class       x         n
## 1      [0,50)   18.84 34.015929
## 2    [50,100)   65.89  3.505522
## 3   [100,500)  158.45  0.785294
## 4 [500,1e+03)  668.27  0.018466
## 5 [1e+03,Inf) 2485.25  0.007766

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.

2.4 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:

employment |> freq_table(age, weights = weights, total = TRUE, print = TRUE)
##        age    x      n
## 1  [15,30) 22.5  26726
## 2  [30,40) 35.0  21045
## 3  [40,50) 45.0  21610
## 4  [50,60) 55.0  20546
## 5 [60,Inf) 65.0  36654
## 6    Total   NA 126582

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

2.5 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.

Table 1: Functions for descriptive statistics
R descstat
mean mean
median median
quantile quantile
var variance
sd stdev
mad madev
modval
medial
gini
skewness
kurtosis

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

z <- wages |> freq_table(wage)
z |> mean()
## [1] 24.12
z |> median()
## [1] 21.94
z |> modval()
##       wage  x  n
## 17 [24,26) 25 69

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,1, we get:

z |> stdev()
## [1] 14.68
z |> variance()
## [1] 215.4
z |> madev()
## [1] 11.57

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"):

z |> quantile(probs = c(0.25, 0.5, 0.75))
## [1] 13.45 21.94 32.31
z |> quantile(y = "mass", probs = c(0.25, 0.5, 0.75))
## [1] 21.58 30.39 46.00

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

z |> median()
## [1] 21.94
z |> medial()
## [1] 30.39

gini computes the Gini coefficient of the series:

z |> gini()
## [1] 0.3403

skewness and kurtosis compute Fisher’s shape statistics:

z |> skewness()
## [1] 1.898
z |> kurtosis()
## [1] 1.48

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

wages |> freq_table(sector) |> modval()
##     sector   n
## 1 services 362

3 Ploting a frequency table

3.1 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),

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

Figure 1: Bar plot for a frequency table

or using a raw table: 

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).

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

Figure 2: Bar plot for integer series

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).

tinyplot(~ sector, wages, type = type_pie())
Figure 3: A pie chart

Figure 3: A pie chart

or:

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).2

tinyplot(~ sector, wages,
         type = type_pie(position = c(1.1, 0.7), init.angle = 10, f = "p",
                         radius = 1.3, pal = "Tableau 10"))
Figure 4: A pie chart with customized initial angle and labels positions

Figure 4: A pie chart with customized initial angle and labels positions

A donut chart is obtained by setting the hole argument to a value strictly between 0 and 1 (see figure 5).

tinyplot(~ sector, wages,
         type = type_pie(position = 1.1, init.angle = 10,
                         radius = 1.3, hole = 0.3, pal = "Set3"))
Figure 5: Donut chart

Figure 5: Donut chart

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.

tinyplot(~ children, data = rgp, type = type_cumul(expand = 0.1, max = 3))
Figure 6: Cummulative distribution

Figure 6: Cummulative distribution

Using a frequency table, we get:

rc <- freq_table(rgp, children, f = "F", max = 3)
tinyplot(F ~ children, data = rc, type = type_cumul(expand = 0.1))

3.2 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),

tinyplot(~ wage, data = wages,
         type = type_histo(breaks = c(10, 20, 30, 40, 50)))
Figure 7: Histogram

Figure 7: Histogram

and descstat::type_freqpoly (see figure 8).

tinyplot(~ wage, data = wages,
         type = type_freqpoly(breaks = c(10, 20, 30, 40, 50)))
Figure 8: Frequency polygons

Figure 8: Frequency polygons

Using a frequency table, the same plots can be obtained using:

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).

tinyplot(~ wage, data = wages, type = type_cumul(expand = 0.05))
Figure 9: Cummulative distribution

Figure 9: Cummulative distribution

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).

tinyplot(~ wage, wages, type = type_lorenz(), col = "lightgrey")
Figure 10: Lorenz curve

Figure 10: Lorenz curve

or, using a frequency table:

ww <- freq_table(wages, wage, f = "FM")
tinyplot(M~ F, ww, type = type_lorenz(), col = "lightgrey")

4 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:

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:

with(wages2, table(wage, size))
##           size
## wage       [0,20) [20,50) [50,100) [100,Inf)
##   [0,10)       74      17       30        63
##   [10,30)     169      75       55       236
##   [30,50)      41      26       21       103
##   [50,Inf)     13       6        7        64

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

with(wages2, table(wage, size)) |> as.data.frame() |> head()
##       wage    size Freq
## 1   [0,10)  [0,20)   74
## 2  [10,30)  [0,20)  169
## 3  [30,50)  [0,20)   41
## 4 [50,Inf)  [0,20)   13
## 5   [0,10) [20,50)   17
## 6  [10,30) [20,50)   75

4.1 Creating contingency tables using cont_table

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

wages2 |> cont_table(wage, size) |> head()
##       wage    size   n
## 1   [0,10)  [0,20)  74
## 2  [10,30)  [0,20) 169
## 3  [30,50)  [0,20)  41
## 4 [50,Inf)  [0,20)  13
## 5   [0,10) [20,50)  17
## 6  [10,30) [20,50)  75

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

wages2 |> cont_table(wage, size) |> wide()
##   wage|size [0,20) [20,50) [50,100) [100,Inf)
## 1    [0,10)     74      17       30        63
## 2   [10,30)    169      75       55       236
## 3   [30,50)     41      26       21       103
## 4  [50,Inf)     13       6        7        64

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.

4.2 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.

tinyplot(wage ~ size, data = wages2,
         type = type_cont.table(size = c(0.5, 5)),
         pch = 16)
Figure 11: Contingency table plot

Figure 11: Contingency table plot

4.3 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.

wht <- wages2 |> cont_table(size, wage)
wht |> joint() |> wide()
##   size|wage [0,10) [10,30) [30,50) [50,Inf)
## 1    [0,20)  0.074   0.169   0.041    0.013
## 2   [20,50)  0.017   0.075   0.026    0.006
## 3  [50,100)  0.030   0.055   0.021    0.007
## 4 [100,Inf)  0.063   0.236   0.103    0.064
wht |> marginal(size)
##        size   x     f
## 1    [0,20)  10 0.297
## 2   [20,50)  35 0.124
## 3  [50,100)  75 0.113
## 4 [100,Inf) 125 0.466
wht |> conditional(size) |> wide()
##   size|wage  [0,10) [10,30) [30,50) [50,Inf)
## 1    [0,20) 0.40217  0.3159  0.2147  0.14444
## 2   [20,50) 0.09239  0.1402  0.1361  0.06667
## 3  [50,100) 0.16304  0.1028  0.1099  0.07778
## 4 [100,Inf) 0.34239  0.4411  0.5393  0.71111

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

4.4 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.

wht |> joint() |> mean()
##    size  wage
## 1 74.03 24.66
wht |> joint() |> stdev()
##    size  wage
## 1 51.16 15.49
wht |> joint() |> variance()
##   size  wage
## 1 2617 240.1
wht |> joint() |> modval()
##   series   value  x     f
## 1   size  [0,20) 10 0.297
## 2   wage [10,30) 20 0.535

The same (univariate) statistics can be obtained using the marginal distribution:

wht |> marginal(size) |> mean()
## [1] 74.03

or even more simply considering the univariate distribution computed by freq_table:

wages2 |> freq_table(size) |> mean()
## [1] 74.03

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:

wht |> conditional(wage) |> mean()
##        size  mean
## 1    [0,20) 20.77
## 2   [20,50) 24.07
## 3  [50,100) 22.21
## 4 [100,Inf) 27.89
wht |> conditional(wage) |> variance()
##        size variance
## 1    [0,20)    180.7
## 2   [20,50)    175.6
## 3  [50,100)    228.3
## 4 [100,Inf)    276.4

4.5 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:

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)
##     om     ev     rv     tv 
##  24.66  10.05 230.03 240.08

Or more simply using the anova method for cont_table objects:

wht_wage <- wht |> anova(wage)
wht_wage
##        size   x     f  mean variance
## 1    [0,20)  10 0.297 20.77    180.7
## 2   [20,50)  35 0.124 24.07    175.6
## 3  [50,100)  75 0.113 22.21    228.3
## 4 [100,Inf) 125 0.466 27.89    276.4

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

wht_wage |> summary()
##   inter intra total   ratio
## 1 10.05   230 240.1 0.04188

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).

tinyplot(wage ~ size, data = wages2, type = type_anova())
Figure 12: Anova plot

Figure 12: Anova plot

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).

4.6 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.

wht |> joint() |> covariance()
## [1] 68.84
wht |> joint() |> correlation()
## [1] 0.08685

The regression line can be computed using the regline function:

rl <- regline(wage ~ size, wht)
rl
## [1] 22.7126  0.0263

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).

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")
Figure 13: Regression line

Figure 13: Regression line

  1. madev computes the mean absolute variation.↩︎

  2. These arguments are the same as graphics::pie.↩︎