Formula versus tidyselect interface

tableone is implemented to allow a subset of columns in a large dataset to be pulled into a table without any fuss. It is also designed with a workflow in mind that involves building statistical models from the data later. We assume the data follows a general pattern in that there are one observation per row, individual columns are specific data points in those observations and may be one of:

• outcome: something that we will be assessing in a statistical model, maybe a continuous outcome, or a time measure, or a logical measure.
• intervention: the thing that is varied between the different observations
• covariates: the other factors that may influence the outcome that we want to control for.

In the end we will want to construct a model that takes the following high level structure:

outcome ~ intervention + covariate_1 + covariate_2 + ... + covariate_n

cgd = survival::cgd %>% 
  # filter to include only the first visit
  dplyr::filter(enum==1) %>% 
  # make the steroids and propylac columns into a logical value
  # see later for a better way of doing this.
  dplyr::mutate(
    steroids = as.logical(steroids),
    propylac = as.logical(propylac)
  )
  

# A basic unstratified population description table is as follows:
cgd %>% describe_population(tidyselect::everything())
#> Warning: Unknown or uninitialised column: `level`.

# define a formula - this might be reused in model building later
formula = Surv(tstart, tstop, status) ~ treat + 
  sex + age + height + weight + inherit + steroids + hos.cat

# set a table relabelling function
rename_cols = function(col) {
  dplyr::case_when(
    col == "hos.cat" ~ "Location",
    col == "steroids" ~ "Steroid treatment",
    TRUE ~ stringr::str_to_sentence(col)
  )
}
options("tableone.labeller"=rename_cols)

# create a simple description
cgd %>% describe_population(formula)

# same as above
formula = Surv(tstart, tstop, status) ~ treat + 
  sex + age + height + weight + inherit + steroids + hos.cat

# labelling function is still active
cgd %>% compare_population(formula)

cgd %>% dplyr::group_by(treat) %>% 
  compare_population(sex,age,height,weight,inherit,steroids,hos.cat)

Analysis of missing data

We need to make sure that not only is the data equivalent between the intervention groups but also that missing data is not unevenly distributed or excessive. Reporting on the frequency of missing data stratified by intervention is also easy, and to demonstrate this we make a data set with 10% of the placebo arm having missing values, but 25% of the treatment arm:

# generate a dataset with values missing not at random compared to the intervention:
cgd_treat = cgd %>% dplyr::mutate(treat = as.character(treat)) %>% dplyr::filter(treat != "placebo")
cgd_placebo = cgd %>% dplyr::mutate(treat = as.character(treat)) %>% dplyr::filter(treat == "placebo")

set.seed(100)
mnar_cgd = dplyr::bind_rows(
  cgd_placebo %>% .make_missing(p_missing = 0.1),
  cgd_treat %>% .make_missing(p_missing = 0.25)
)

Comparing this new data set we see that there is significant differences in some of the data (but not the steroids variable). As this is quite a small dataset it is not sufficiently powered to reliably detect the difference in missingness at this level (15% difference).

# compare the MNAR dataset against the intervention:
formula = Surv(tstart, tstop, status) ~ treat + 
  sex + age + height + weight + inherit + steroids + hos.cat

mnar_cgd %>% compare_missing(formula)

with this analysis it is useful to be able to update the analysis formula removing the variables with missing data so that we are confident the models are based on reasonable data.

# formula can also be a list of formulae
new_formula = mnar_cgd %>% tableone::remove_missing(formula)
#> More than 10% of data is missing for variables Sex, Age, Height, Weight, Inherit, Steroid treatment, Location.
#> Data is missing not at random (compared to Treat) at a p-value<0.007 (0.05 over 7 comparisons) for variables Age.

print(new_formula)
#> [[1]]
#> Surv(tstart, tstop, status) ~ treat

Conversion of discrete data

Using this new data set with missing data it may be necessary to discretise some or all of the data, or convert logical values into properly named factors.

decade = function(x) sprintf("%d-%d",x-(x%%10),x-(x%%10)+9)

discrete_cgd = mnar_cgd %>% 
  # pick out the first episode
  dplyr::filter(enum == 1) %>%
  # convert data
  make_factors(
    steroids,propylac,age,weight,height,
    .logical = c("received","not received"),
    .numeric = list(
      age="{decade(value)}",
      weight="{ifelse(value<20,'<20','20+')}",
      height="{ifelse(value<mean(value, na.rm=TRUE),'below average','above average')}"
    )
  )

formula = Surv(tstart, tstop, status) ~ treat + 
  sex + age + height + weight + inherit + steroids + hos.cat


old = options("tableone.show_pvalue_method"=TRUE)
# This comparison implicitly ignores missing values.
t = discrete_cgd %>% compare_population(formula)
options(old)

t

# N.B. The following option is involved when converting integer data
# which decides how many levels of integer data are considered discrete
# and when to decide integer data can be treated as continuous:
options("tableone.max_discrete_levels"=0)
# and is described in the documentation for make_factors().

Making missing factors explicit:

In the comparison above missing values were not included, and we should be cautious of the findings. Because of the missingness tableone will not calculate p-values. If factor values are missing (as in this case) then we can include them as a new group and get a more robust comparison which includes the distribution of missingness, and for which we can calculate a p-value. However previously ordered variables, are now regarded as unordered as we cannot determine the value of a missing level.

discrete_cgd %>% explicit_na() %>% compare_population(formula)

Non biomedical data

Beyond the bio-medical example tableone can make any more general comparison between data that has a structure like:

~ group + observation_1 + observation_2 + ... + observation_n

We will use the iris and the diamonds datasets to demonstrate this more general use case for tableone.

# revert the labeller setting to the default
# and additionally hide the footer.
old = options(
  "tableone.labeller"=NULL,
  "tableone.show_pvalue_method"=FALSE,
  "tableone.hide_footer"=TRUE)

# the heuristics detect that Petals in the iris data set are not normally
# distributed and hence report median and IQR:
iris %>% dplyr::group_by(Species) %>% compare_population(tidyselect::everything())

options(old)

The missing_diamonds data set which is included in this package has 10% of the values removed. This demonstrates the need for reporting the denominator.

# The counts sometimes seem redundant if there is no missing information:
# however in a data set with missing values the denominators are important:
missing_diamonds %>% describe_population(tidyselect::everything())