2 Descriptive statistics
Many options currently exists for displaying descriptive statistics. I recommend the function summarise from the papeR package. It does a great job for both continuous and categorical (factors) outcomes.
The estat function from the pubh package displays descriptive statistics of continuous outcomes and like summarise, it can use labels to display nice tables. estat has less options than summarise and as such, is simpler to use, which is the main intention of the pubh package. Furthermore, 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 variance between groups (factors).
Some examples. We start by loading packages as suggested in the Template of this package.
library(effects, warn.conflicts = FALSE)
library(latex2exp, warn.conflicts = FALSE)
library(multcomp, warn.conflicts = FALSE)
library(pander, warn.conflicts = FALSE)
library(papeR, warn.conflicts = FALSE)
library(pubh, warn.conflicts = FALSE)
panderOptions("table.split.table", Inf)
set.alignment("right", row.names = "left", permanent = TRUE)
trellis.par.set(tactile.theme())We will use a data set about blood counts of T cells from patients in remission from Hodgkin’s disease or in remission from disseminated malignancies.
data(Hodgkin)
Hodgkin$Ratio <- Hodgkin$CD4/Hodgkin$CD8
pander(head(Hodgkin))| CD4 | CD8 | Group | Ratio |
|---|---|---|---|
| 396 | 836 | Hodgkin | 0.4737 |
| 568 | 978 | Hodgkin | 0.5808 |
| 1212 | 1678 | Hodgkin | 0.7223 |
| 171 | 212 | Hodgkin | 0.8066 |
| 554 | 670 | Hodgkin | 0.8269 |
| 1104 | 1335 | Hodgkin | 0.827 |
Descriptive statistics for CD4+ T cells:
estat(~ CD4, data = Hodgkin)
#> N Min. Max. Mean Median SD CV
#> 1 CD4 40 116 2415 672.62 528.5 470.49 0.7In the previous code, the left-hand side of the formula is empty as it’s the case when working with only one variable.
We can use pander to display nice tables with labels.
pander(estat(~ CD4, data = Hodgkin, label = "CD4^+^ T cells (cells / mm^3^)"))| N | Min. | Max. | Mean | Median | SD | CV | |
|---|---|---|---|---|---|---|---|
| CD4+ T cells (cells / mm3) | 40 | 116 | 2415 | 672.6 | 528.5 | 470.5 | 0.7 |
In the previous command, we used the Markup language for the superscripts in the label. For subscripts, the syntax is (for example): H~2~O.
For stratification, estat recognises the following two syntaxes:
outcome ~ categorical predictor~ outcome | categorical predictor
For example:
pander(estat(~ Ratio|Group, data = Hodgkin, label = "CD4^+^/CD8^+^ T cells"))| Group | N | Min. | Max. | Mean | Median | SD | CV | |
|---|---|---|---|---|---|---|---|---|
| CD4+/CD8+ T cells | Non-Hodgkin | 20 | 1.1 | 3.49 | 2.12 | 2.15 | 0.73 | 0.34 |
| Hodgkin | 20 | 0.47 | 3.82 | 1.5 | 1.19 | 0.91 | 0.61 |
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.6293991What about normality? We can look at the QQ-plots against the standard Normal distribution.
lab <- TeX("CD4$^+$ / CD8$^+$ T cells")
qq_plot(~ Ratio|Group, data = Hodgkin, ylab = lab)
Let’s say we choose a non-parametric test to compare the groups:
wilcox.test(Ratio ~ Group, data = Hodgkin)
#>
#> Wilcoxon rank sum test
#>
#> data: Ratio by Group
#> W = 298, p-value = 0.007331
#> alternative hypothesis: true location shift is not equal to 02.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:
strip_error(Ratio ~ Group, data = Hodgkin, ylab = lab)
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. With the latticeExtra package (loaded by default with pubh) allows the use of layers. We will draw a line using panel.segments and text with panel.text.
For panel.segments we provide the coordinates of the two points at the ends of the line in the following order: x1, y1, x2, y2. x coordinates are easy: from 1 to 2. For the y coordinates, we choose a value higher than the maximum observation (from our descriptives, higher than 3.82). We also need to adjust the limits of the y-axis to accommodate for our line and text. It could be some trial and error one has the plot that is more visually appealing.
The commands are as follows:
strip_error(Ratio ~ Group, data = Hodgkin, ylab = lab, ylim = c(0, 4.5)) +
layer(panel.segments(1, 4, 2, 4, lwd=1.5)) +
layer(panel.text(1.5, 4.1, "**"))
For larger samples we would use bar charts with error bars. For example:
data(birthwt)
birthwt$smoke <- factor(birthwt$smoke, labels = c("Non-smoker", "Smoker"))
birthwt$race <- factor(birthwt$race, labels = c("White", "African American", "Other"))Bar charts should always start at the origin, the upper limit on the \(y\) scale should be higher than the maximum upper confidence limit. For example:
pander(gen_bst_df(bwt ~ smoke, data = birthwt))| bwt | LowerCI | UpperCI | smoke |
|---|---|---|---|
| 3056 | 2923 | 3186 | Non-smoker |
| 2772 | 2616 | 2910 | Smoker |
bar_error(bwt ~ smoke, data = birthwt, ylab = "Birth weight (g)", ylim = c(0, 3500))
Quick normality check:
qq_plot(~ bwt|smoke, data = birthwt, ylab = "Birth weight (g)")
The (unadjusted) 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.919The final plot showing (unadjusted) difference:
bar_error(bwt ~ smoke, data = birthwt, ylab = "Birth weight (g)", ylim = c(0, 3500)) +
layer(panel.segments(1, 3300, 2, 3300, lwd=1.5)) +
layer(panel.text(1.5, 3350, "**"))
Both strip_error and bar_error can generate plots stratified by a third variable, for example:
bar_error(bwt ~ smoke|race, data = birthwt, ylab = "Birth weight (g)", ylim = c(0, 3800),
col = c("gray80", "gray30"))
strip_error(bwt ~ smoke|race, data = birthwt, ylab = "Birth weight (g)", cex = 0.3)
