2 Descriptive statistics
Many options currently exists for displaying descriptive statistics. I recommend the function mytable from the moonBook package for constructing tables of descriptive statistics where results don’t need to be stratified.
In Public Health and Epidemiology it is common to have a categorical outcome and thus, report descriptive statistics stratified by the outcome. Package pubh introduces the function cross_tab as a wrapper to functions from finalfit. The idea is to construct these tables, using formulas, in a simple way.
The estat function from the pubh package displays descriptive statistics of continuous outcomes and like mytable, it can use labels to display nice tables. estat. As a way to aid in the understanding of the variability, estat displays also the relative dispersion (coefficient of variation) which is of particular interest for comparing variances between groups (factors).
Some examples. We start by loading packages.
rm(list = ls())
library(car)
library(broom)
library(tidyverse)
library(ggfortify)
library(mosaic)
library(huxtable)
library(jtools)
library(latex2exp)
library(pubh)
library(sjlabelled)
library(sjPlot)
library(sjmisc)
theme_set(sjPlot::theme_sjplot2(base_size = 10))
theme_update(legend.position = "top")
# options('huxtable.knit_print_df' = FALSE)
options('huxtable.knit_print_df_theme' = theme_article)
options('huxtable.autoformat_number_format' = list(numeric = "%5.2f"))
knitr::opts_chunk$set(collapse = TRUE, comment = NA)We will use a data set about a study of onchocerciasis in Sierra Leone.
| id | mf | area | agegrp | sex | mfload | lesions |
|---|---|---|---|---|---|---|
| 1 | Infected | Savannah | 20-39 | Female | 1 | No |
| 2 | Infected | Rainforest | 40+ | Male | 3 | No |
| 3 | Infected | Savannah | 40+ | Female | 1 | No |
| 4 | Not-infected | Rainforest | 20-39 | Female | 0 | No |
| 5 | Not-infected | Savannah | 40+ | Female | 0 | No |
| 6 | Not-infected | Rainforest | 20-39 | Female | 0 | No |
A two-by-two contingency table:
Oncho %>%
mutate(
mf = relevel(mf, ref = "Infected")
) %>%
copy_labels(Oncho) %>%
cross_tab(mf ~ area) %>%
theme_article()| Infection | |||
|---|---|---|---|
| Infected | Not-infected | Total | |
| (N=822) | (N=480) | (N=1302) | |
| Residence | |||
| - Savannah | 281 (34.2%) | 267 (55.6%) | 548 (42.1%) |
| - Rainforest | 541 (65.8%) | 213 (44.4%) | 754 (57.9%) |
Table with all descriptive statistics except id and mfload:
Oncho %>%
select(- c(id, mfload)) %>%
mutate(
mf = relevel(mf, ref = "Infected")
) %>%
copy_labels(Oncho) %>%
cross_tab(mf ~ area +.) %>%
theme_article()| Infection | |||
|---|---|---|---|
| Infected | Not-infected | Total | |
| (N=822) | (N=480) | (N=1302) | |
| Residence | |||
| - Savannah | 281 (34.2%) | 267 (55.6%) | 548 (42.1%) |
| - Rainforest | 541 (65.8%) | 213 (44.4%) | 754 (57.9%) |
| Age group (years) | |||
| - 5-9 | 46 ( 5.6%) | 156 (32.5%) | 202 (15.5%) |
| - 10-19 | 99 (12.0%) | 119 (24.8%) | 218 (16.7%) |
| - 20-39 | 299 (36.4%) | 125 (26.0%) | 424 (32.6%) |
| - 40+ | 378 (46.0%) | 80 (16.7%) | 458 (35.2%) |
| Sex | |||
| - Male | 426 (51.8%) | 190 (39.6%) | 616 (47.3%) |
| - Female | 396 (48.2%) | 290 (60.4%) | 686 (52.7%) |
| Severe eye lesions? | |||
| - No | 640 (77.9%) | 461 (96.0%) | 1101 (84.6%) |
| - Yes | 182 (22.1%) | 19 ( 4.0%) | 201 (15.4%) |
Next, we use a data set about blood counts of T cells from patients in remission from Hodgkin’s disease or in remission from disseminated malignancies. We generate the new variable Ratio and add labels to the variables.
data(Hodgkin)
Hodgkin <- Hodgkin %>%
mutate(Ratio = CD4/CD8) %>%
var_labels(
Ratio = "CD4+ / CD8+ T-cells ratio"
)
Hodgkin %>% head()| CD4 | CD8 | Group | Ratio |
|---|---|---|---|
| 396 | 836 | Hodgkin | 0.47 |
| 568 | 978 | Hodgkin | 0.58 |
| 1212 | 1678 | Hodgkin | 0.72 |
| 171 | 212 | Hodgkin | 0.81 |
| 554 | 670 | Hodgkin | 0.83 |
| 1104 | 1335 | Hodgkin | 0.83 |
Descriptive statistics for CD4+ T cells:
| N | Min. | Max. | Mean | Median | SD | CV | |
|---|---|---|---|---|---|---|---|
| CD4+ T-cells | 40.00 | 116.00 | 2415.00 | 672.62 | 528.50 | 470.49 | 0.70 |
In the previous code, the left-hand side of the formula is empty as it’s the case when working with only one variable.
For stratification, estat recognises the following two syntaxes:
outcome ~ exposure~ outcome | exposure
where, outcome is continuous and exposure is a categorical (factor) variable.
For example:
| Disease | N | Min. | Max. | Mean | Median | SD | CV | |
|---|---|---|---|---|---|---|---|---|
| CD4+ / CD8+ T-cells ratio | Non-Hodgkin | 20.00 | 1.10 | 3.49 | 2.12 | 2.15 | 0.73 | 0.34 |
| Hodgkin | 20.00 | 0.47 | 3.82 | 1.50 | 1.19 | 0.91 | 0.61 |
As before, we can report a table of descriptive statistics for all variables in the data set:
Hodgkin %>%
mutate(
Group = relevel(Group, ref = "Hodgkin")
) %>%
copy_labels(Hodgkin) %>%
cross_tab(Group ~ CD4 + ., method = 2) %>%
theme_article() %>%
add_footnote("Values are medians with interquartile range.")| Disease | |||
|---|---|---|---|
| Hodgkin | Non-Hodgkin | Total | |
| (N=20) | (N=20) | (N=40) | |
| CD4+ T-cells | 681.5 [396.5;1158.0] | 433.0 [345.0;718.0] | 528.5 [375.0;930.0] |
| CD8+ T-cells | 447.5 [298.5;823.5] | 231.5 [146.5;325.0] | 319.0 [206.0;601.0] |
| CD4+ / CD8+ T-cells ratio | 1.2 [ 0.8; 2.0] | 2.2 [ 1.6; 2.7] | 1.7 [ 1.1; 2.4] |
| Values are medians with interquartile range. | |||
2.1 Inferential statistics
From the descriptive statistics of Ratio we know that the relative dispersion in the Hodgkin group is almost as double as the relative dispersion in the Non-Hodgkin group.
For new users of R it helps to have a common syntax in most of the commands, as R could be challenging and even intimidating at times. We can test if the difference in variance is statistically significant with the var.test command, which uses can use a formula syntax:
var.test(Ratio ~ Group, data = Hodgkin)
F test to compare two variances
data: Ratio by Group
F = 0.6294, num df = 19, denom df = 19, p-value = 0.3214
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.2491238 1.5901459
sample estimates:
ratio of variances
0.6293991 What about normality? We can look at the QQ-plots against the standard Normal distribution.
Let’s say we choose a non-parametric test to compare the groups:
2.2 Graphical output
For relatively small samples (for example, less than 30 observations per group) is a standard practice to show the actual data in dot plots with error bars. The pubh package offers two options to show graphically differences in continuous outcomes among groups:
- For small samples:
strip_error - For medium to large samples:
bar_error
For our current example:
In the previous code, axis_labs applies labels from labelled data to the axis.
The error bars represent 95% confidence intervals around mean values.
Is relatively easy to add a line on top, to show that the two groups are significantly different. The function gf_star needs the reference point on how to draw an horizontal line to display statistical difference or to annotate a plot; in summary, gf_star:
- Draws an horizontal line from \((x1, y2)\) to \((x2, y2)\).
- Draws a vertical line from \((x1, y1)\) to \((x1, y1)\).
- Draws a vertical line from \((x2, y1)\) to \((x2, y1)\).
- Writes a character (the legend or “star”) at the mid point between \(x1\) and \(x2\) at high \(y3\).
Thus:
\[ y1 < y2 < y3 \]
In our current example:
Hodgkin %>%
strip_error(Ratio ~ Group) %>%
axis_labs() %>%
gf_star(x1 = 1, y1 = 4, x2 = 2, y2 = 4.05, y3 = 4.1, "**")
For larger samples we could use bar charts with error bars. For example:
data(birthwt, package = "MASS")
birthwt <- birthwt %>%
mutate(
smoke = factor(smoke, labels = c("Non-smoker", "Smoker")),
Race = factor(race > 1, labels = c("White", "Non-white")),
race = factor(race, labels = c("White", "Afican American", "Other"))
) %>%
var_labels(
bwt = 'Birth weight (g)',
smoke = 'Smoking status',
race = 'Race',
)
Quick normality check:
The (unadjusted) \(t\)-test:
t.test(bwt ~ smoke, data = birthwt)
Welch Two Sample t-test
data: bwt by smoke
t = 2.7299, df = 170.1, p-value = 0.007003
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
78.57486 488.97860
sample estimates:
mean in group Non-smoker mean in group Smoker
3055.696 2771.919 The final plot with annotation to highlight statistical difference (unadjusted):
birthwt %>%
bar_error(bwt ~ smoke) %>%
axis_labs() %>%
gf_star(x1 = 1, x2 = 2, y1 = 3400, y2 = 3500, y3 = 3550, "**")
Both strip_error and bar_error can generate plots stratified by a third variable, for example:
birthwt %>%
bar_error(bwt ~ smoke, fill = ~ Race) %>%
gf_refine(ggsci::scale_fill_jama()) %>%
axis_labs()
birthwt %>%
strip_error(bwt ~ smoke, pch = ~ Race, col = ~ Race) %>%
gf_refine(ggsci::scale_color_jama()) %>%
axis_labs()
