All tables examples

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

All tables examples

Ewen Harrison

1 Cross tables

Two-way tables are used extensively in healthcare research, e.g. a 2x2 table comparing two factors with two levels each, or table 1 from a typical clinical study or trial

The main functions all take a dependent variable - the outcome - and explanatory variables - predictors or exposures (any number categorical or continuous variables).

1.01 Default

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
  summary_factorlist(dependent, explanatory) -> t

label	levels	No	Yes
Age (years)	Mean (SD)	59.8 (11.9)	58.4 (13.3)
Age	<40 years	68 (7.5)	2 (7.4)
	40-59 years	334 (37.0)	10 (37.0)
	60+ years	500 (55.4)	15 (55.6)
Sex	Female	432 (47.9)	13 (48.1)
	Male	470 (52.1)	14 (51.9)
Obstruction	No	715 (81.2)	17 (63.0)
	Yes	166 (18.8)	10 (37.0)

Note, chi-squared warnings will be generated when the expected count in any cell is less than 5. Fisher’s exact test can be used as below, or go straight to a univariable logistic regression, e.g. colon_s %>% finalfit(dependent, explanatory)

1.02 Add or edit variable labels

library(finalfit)
library(dplyr)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    mutate(
        sex.factor = ff_label(sex.factor, "Gender")
    ) %>% 
  summary_factorlist(dependent, explanatory) -> t

label	levels	No	Yes
Age (years)	Mean (SD)	59.8 (11.9)	58.4 (13.3)
Age	<40 years	68 (7.5)	2 (7.4)
	40-59 years	334 (37.0)	10 (37.0)
	60+ years	500 (55.4)	15 (55.6)
Gender	Female	432 (47.9)	13 (48.1)
	Male	470 (52.1)	14 (51.9)
Obstruction	No	715 (81.2)	17 (63.0)
	Yes	166 (18.8)	10 (37.0)

1.03 P-value for hypothesis test

Defaults are chi-squared for categorical explanatory variables and an F-test for continuous (aov, analysis of variance). Alternatives can be specified as per below.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE) -> t

label	levels	No	Yes	p
Age (years)	Mean (SD)	59.8 (11.9)	58.4 (13.3)	0.542
Age	<40 years	68 (7.5)	2 (7.4)	1.000
	40-59 years	334 (37.0)	10 (37.0)
	60+ years	500 (55.4)	15 (55.6)
Sex	Female	432 (47.9)	13 (48.1)	1.000
	Male	470 (52.1)	14 (51.9)
Obstruction	No	715 (81.2)	17 (63.0)	0.035
	Yes	166 (18.8)	10 (37.0)

1.04 With Fisher’s exact test

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
  summary_factorlist(dependent, explanatory, p = TRUE, p_cat = "fisher") -> t

label	levels	No	Yes	p
Age (years)	Mean (SD)	59.8 (11.9)	58.4 (13.3)	0.542
Age	<40 years	68 (7.5)	2 (7.4)	1.000
	40-59 years	334 (37.0)	10 (37.0)
	60+ years	500 (55.4)	15 (55.6)
Sex	Female	432 (47.9)	13 (48.1)	1.000
	Male	470 (52.1)	14 (51.9)
Obstruction	No	715 (81.2)	17 (63.0)	0.026
	Yes	166 (18.8)	10 (37.0)

1.05 Parametric explanatory variables

Summaries for continuous explanatory variables are mean (standard deviation) with aov statistical test by default. The statistical test can be changed to the Welch t-test when there are two dependent variable levels if desired.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
  summary_factorlist(dependent, explanatory, p = TRUE, p_cont_para = "t.test") -> t
#> Warning: There was 1 warning in `dplyr::summarise()`.
#> ℹ In argument: `chisq.test(age.factor, perfor.factor)$p.value`.
#> Caused by warning in `chisq.test()`:
#> ! Chi-squared approximation may be incorrect

label	levels	No	Yes	p
Age (years)	Mean (SD)	59.8 (11.9)	58.4 (13.3)	0.586
Age	<40 years	68 (7.5)	2 (7.4)	1.000
	40-59 years	334 (37.0)	10 (37.0)
	60+ years	500 (55.4)	15 (55.6)
Sex	Female	432 (47.9)	13 (48.1)	1.000
	Male	470 (52.1)	14 (51.9)
Obstruction	No	715 (81.2)	17 (63.0)	0.035
	Yes	166 (18.8)	10 (37.0)

1.06 Non-parametric explanatory variables

If desired, all continuous explanatory variables can be considered non-parametric. Summaries will be median (interquartile range) and the statistical test is Kruskal-Wallis/Mann-Whitney U. Use cont_range = FALSE if wish single-digit IQR, i.e. Q3-Q1.

library(finalfit)
explanatory = c("age", "nodes", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, cont = "median") -> t

label	levels	No	Yes	p
Age (years)	Median (IQR)	61.0 (53.0 to 69.0)	60.0 (50.0 to 68.0)	0.578
nodes	Median (IQR)	2.0 (1.0 to 5.0)	3.0 (2.0 to 4.0)	0.125
Age	<40 years	68 (7.5)	2 (7.4)	1.000
	40-59 years	334 (37.0)	10 (37.0)
	60+ years	500 (55.4)	15 (55.6)
Sex	Female	432 (47.9)	13 (48.1)	1.000
	Male	470 (52.1)	14 (51.9)
Obstruction	No	715 (81.2)	17 (63.0)	0.035
	Yes	166 (18.8)	10 (37.0)

1.07 Select specific non-parametric variables

Many have asked in the past if only particular variables can be considered non-parametric. The argument cont_nonpara can take a vector (e.g. c(1, 2, 3, 4)) of values corresponding to the explanatory variable to specify which should be summarised as a median and be passed to a non-parametric test.

library(finalfit)
explanatory = c("age", "nodes", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, cont_nonpara = c(2)) -> t

label	levels	No	Yes	p
Age (years)	Mean (SD)	59.8 (11.9)	58.4 (13.3)	0.542
nodes	Median (IQR)	2.0 (1.0 to 5.0)	3.0 (2.0 to 4.0)	0.125
Age	<40 years	68 (7.5)	2 (7.4)	1.000
	40-59 years	334 (37.0)	10 (37.0)
	60+ years	500 (55.4)	15 (55.6)
Sex	Female	432 (47.9)	13 (48.1)	1.000
	Male	470 (52.1)	14 (51.9)
Obstruction	No	715 (81.2)	17 (63.0)	0.035
	Yes	166 (18.8)	10 (37.0)

1.08 Missing values for the explanatory variables

Always consider summarising missing values when describing your data.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, na_include = TRUE) -> t

label	levels	No	Yes	p
Age (years)	Mean (SD)	59.8 (11.9)	58.4 (13.3)	0.542
Age	<40 years	68 (7.5)	2 (7.4)	1.000
	40-59 years	334 (37.0)	10 (37.0)
	60+ years	500 (55.4)	15 (55.6)
	(Missing)	0 (0.0)	0 (0.0)
Sex	Female	432 (47.9)	13 (48.1)	1.000
	Male	470 (52.1)	14 (51.9)
	(Missing)	0 (0.0)	0 (0.0)
Obstruction	No	715 (79.3)	17 (63.0)	0.035
	Yes	166 (18.4)	10 (37.0)
	(Missing)	21 (2.3)	0 (0.0)

1.09 Pass missing values to statistical tests

This is a change from the default behaviour introduced in Finalfit 1.0.0. Previously, when missing data was presented it was also considered as a level in the statistical test. This may or may not be desired. Control this now using na_to_p = TRUE to include missing data in test. A message is produced reminding you that you are doing that.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, na_include = TRUE, 
                                         na_to_p = TRUE) -> t
#> Explanatory variable(s) missing data included in hypothesis test (p-value).

label	levels	No	Yes	p
Age (years)	Mean (SD)	59.8 (11.9)	58.4 (13.3)	0.542
Age	<40 years	68 (7.5)	2 (7.4)	1.000
	40-59 years	334 (37.0)	10 (37.0)
	60+ years	500 (55.4)	15 (55.6)
	(Missing)	0 (0.0)	0 (0.0)
Sex	Female	432 (47.9)	13 (48.1)	1.000
	Male	470 (52.1)	14 (51.9)
	(Missing)	0 (0.0)	0 (0.0)
Obstruction	No	715 (79.3)	17 (63.0)	0.042
	Yes	166 (18.4)	10 (37.0)
	(Missing)	21 (2.3)	0 (0.0)

1.10 Row proportions (rather than column)

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, na_include = TRUE,
                                         column = FALSE) -> t

label	levels	No	Yes	p
Age (years)	Mean (SD)	59.8 (11.9)	58.4 (13.3)	0.542
Age	<40 years	68 (97.1)	2 (2.9)	1.000
	40-59 years	334 (97.1)	10 (2.9)
	60+ years	500 (97.1)	15 (2.9)
	(Missing)	0 (NaN)	0 (NaN)
Sex	Female	432 (97.1)	13 (2.9)	1.000
	Male	470 (97.1)	14 (2.9)
	(Missing)	0 (NaN)	0 (NaN)
Obstruction	No	715 (97.7)	17 (2.3)	0.035
	Yes	166 (94.3)	10 (5.7)
	(Missing)	21 (100.0)	0 (0.0)

1.11 Total column

The terms total column was introduced before this function summarised continuous variables. It would be better to be “All data” or something similar, as the continous explanatory variables a summary statistic is produced for all data. However, to keep backwards compatibility we have left it unchanged for now. For producing row totals including continous explanatory variables, see add_row_total below.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, cont = "median", na_include = TRUE,
                                         column = TRUE, total_col = TRUE) -> t

label	levels	No	Yes	Total	p
Age (years)	Median (IQR)	61.0 (53.0 to 69.0)	60.0 (50.0 to 68.0)	61.0 (53.0 to 69.0)	0.578
Age	<40 years	68 (7.5)	2 (7.4)	70 (7.5)	1.000
	40-59 years	334 (37.0)	10 (37.0)	344 (37.0)
	60+ years	500 (55.4)	15 (55.6)	515 (55.4)
	(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	Female	432 (47.9)	13 (48.1)	445 (47.9)	1.000
	Male	470 (52.1)	14 (51.9)	484 (52.1)
	(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	No	715 (79.3)	17 (63.0)	732 (78.8)	0.035
	Yes	166 (18.4)	10 (37.0)	176 (18.9)
	(Missing)	21 (2.3)	0 (0.0)	21 (2.3)

1.12 Row totals with missing

This was introduced to deal with the problem of summarising missing data for continuous variables. By default, it provides the total N for the variable and includes a column enumerating missing values.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, na_include = TRUE,
                                         total_col = TRUE,
                                         add_row_total = TRUE) -> t

label	Total N	Missing N	levels	No	Yes	Total	p
Age (years)	929 (100.0)	0	Mean (SD)	59.8 (11.9)	58.4 (13.3)	59.8 (11.9)	0.542
Age	929 (100.0)	0	<40 years	68 (7.5)	2 (7.4)	70 (7.5)	1.000
			40-59 years	334 (37.0)	10 (37.0)	344 (37.0)
			60+ years	500 (55.4)	15 (55.6)	515 (55.4)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	929 (100.0)	0	Female	432 (47.9)	13 (48.1)	445 (47.9)	1.000
			Male	470 (52.1)	14 (51.9)	484 (52.1)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	908 (97.7)	21	No	715 (79.3)	17 (63.0)	732 (78.8)	0.035
			Yes	166 (18.4)	10 (37.0)	176 (18.9)
			(Missing)	21 (2.3)	0 (0.0)	21 (2.3)

1.13 Row totals without missing

Remove missing column.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, na_include = TRUE,
                                         total_col = TRUE,
                                         add_row_total = TRUE,
                                         include_row_missing_col = FALSE) -> t

label	Total N	levels	No	Yes	Total	p
Age (years)	929 (100.0)	Mean (SD)	59.8 (11.9)	58.4 (13.3)	59.8 (11.9)	0.542
Age	929 (100.0)	<40 years	68 (7.5)	2 (7.4)	70 (7.5)	1.000
		40-59 years	334 (37.0)	10 (37.0)	344 (37.0)
		60+ years	500 (55.4)	15 (55.6)	515 (55.4)
		(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	929 (100.0)	Female	432 (47.9)	13 (48.1)	445 (47.9)	1.000
		Male	470 (52.1)	14 (51.9)	484 (52.1)
		(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	908 (97.7)	No	715 (79.3)	17 (63.0)	732 (78.8)	0.035
		Yes	166 (18.4)	10 (37.0)	176 (18.9)
		(Missing)	21 (2.3)	0 (0.0)	21 (2.3)

1.14 Row totals with user-specified column names

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE,
                                         total_col = TRUE,
                                         add_row_total = TRUE,
                                         row_totals_colname = "N (total)",
                                         row_missing_colname = "N (missing)") -> t

label	N (total)	N (missing)	levels	No	Yes	Total	p
Age (years)	929 (100.0)	0	Mean (SD)	59.8 (11.9)	58.4 (13.3)	59.8 (11.9)	0.542
Age	929 (100.0)	0	<40 years	68 (7.5)	2 (7.4)	70 (7.5)	1.000
			40-59 years	334 (37.0)	10 (37.0)	344 (37.0)
			60+ years	500 (55.4)	15 (55.6)	515 (55.4)
Sex	929 (100.0)	0	Female	432 (47.9)	13 (48.1)	445 (47.9)	1.000
			Male	470 (52.1)	14 (51.9)	484 (52.1)
Obstruction	908 (97.7)	21	No	715 (81.2)	17 (63.0)	732 (80.6)	0.035
			Yes	166 (18.8)	10 (37.0)	176 (19.4)

1.15 Order a variable by total

This is intended for when there is only one explanatory variable.

library(finalfit)
explanatory = c("extent.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, cont = "median", na_include = TRUE,
                                         column = TRUE, total_col = TRUE, orderbytotal = TRUE) -> t

label	levels	No	Yes	Total	p
Extent of spread	Serosa	736 (81.6)	23 (85.2)	759 (81.7)	0.200
	Muscle	105 (11.6)	1 (3.7)	106 (11.4)
	Adjacent structures	40 (4.4)	3 (11.1)	43 (4.6)
	Submucosa	21 (2.3)	0 (0.0)	21 (2.3)
	(Missing)	0 (0.0)	0 (0.0)	0 (0.0)

1.17 Add column totals

Column totals can be added, and by default are presented with a row percentage.

explanatory = c("age.factor", "sex.factor")
dependent = "rx.factor"

colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE,
                                         add_col_totals = TRUE) -> t

label	levels	Obs	Lev	Lev+5FU	p
Total N (%)		315 (33.9)	310 (33.4)	304 (32.7)
Age	<40 years	25 (7.9)	19 (6.1)	26 (8.6)	0.572
	40-59 years	124 (39.4)	115 (37.1)	105 (34.5)
	60+ years	166 (52.7)	176 (56.8)	173 (56.9)
Sex	Female	149 (47.3)	133 (42.9)	163 (53.6)	0.028
	Male	166 (52.7)	177 (57.1)	141 (46.4)

1.18 Add column totals without proportion.

explanatory = c("age.factor", "sex.factor")
dependent = "rx.factor"

colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE,
                                         add_col_totals = TRUE,
                                         include_col_totals_percent = FALSE)    -> t

label	levels	Obs	Lev	Lev+5FU	p
Total N		315	310	304
Age	<40 years	25 (7.9)	19 (6.1)	26 (8.6)	0.572
	40-59 years	124 (39.4)	115 (37.1)	105 (34.5)
	60+ years	166 (52.7)	176 (56.8)	173 (56.9)
Sex	Female	149 (47.3)	133 (42.9)	163 (53.6)	0.028
	Male	166 (52.7)	177 (57.1)	141 (46.4)

1.19 Add column totals with user-specified row name and prefix.

explanatory = c("age.factor", "sex.factor")
dependent = "rx.factor"

colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE,
                                         add_col_totals = TRUE,
                                         include_col_totals_percent = FALSE,
                                         col_totals_rowname = "",
                                         col_totals_prefix = "N=")  -> t

label	levels	Obs	Lev	Lev+5FU	p
		N=315	N=310	N=304
Age	<40 years	25 (7.9)	19 (6.1)	26 (8.6)	0.572
	40-59 years	124 (39.4)	115 (37.1)	105 (34.5)
	60+ years	166 (52.7)	176 (56.8)	173 (56.9)
Sex	Female	149 (47.3)	133 (42.9)	163 (53.6)	0.028
	Male	166 (52.7)	177 (57.1)	141 (46.4)

1.20 Label with `dependent` name

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    summary_factorlist(dependent, explanatory, p = TRUE, cont = "median", na_include = TRUE,
                                         column = TRUE, total_col = TRUE, add_dependent_label = TRUE) -> t

Dependent: Perforation		No	Yes	Total	p
Age (years)	Median (IQR)	61.0 (53.0 to 69.0)	60.0 (50.0 to 68.0)	61.0 (53.0 to 69.0)	0.578
Age	<40 years	68 (7.5)	2 (7.4)	70 (7.5)	1.000
	40-59 years	334 (37.0)	10 (37.0)	344 (37.0)
	60+ years	500 (55.4)	15 (55.6)	515 (55.4)
	(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	Female	432 (47.9)	13 (48.1)	445 (47.9)	1.000
	Male	470 (52.1)	14 (51.9)	484 (52.1)
	(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	No	715 (79.3)	17 (63.0)	732 (78.8)	0.035
	Yes	166 (18.4)	10 (37.0)	176 (18.9)
	(Missing)	21 (2.3)	0 (0.0)	21 (2.3)

The dependent name cannot be passed directly to the table intentionally. This is to avoid errors when code is copied and the name is not updated. Change the dependent label using the following. The prefix (“Dependent:”) and any suffix can be altered.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
  dplyr::mutate(
    perfor.factor = ff_label(perfor.factor, "Perforated cancer")
    ) %>% 
  summary_factorlist(dependent, explanatory, p = TRUE, cont = "median", na_include = TRUE,
    column = TRUE, total_col = TRUE, add_dependent_label = TRUE, dependent_label_prefix = "") -> t

Perforated cancer		No	Yes	Total	p
Age (years)	Median (IQR)	61.0 (53.0 to 69.0)	60.0 (50.0 to 68.0)	61.0 (53.0 to 69.0)	0.578
Age	<40 years	68 (7.5)	2 (7.4)	70 (7.5)	1.000
	40-59 years	334 (37.0)	10 (37.0)	344 (37.0)
	60+ years	500 (55.4)	15 (55.6)	515 (55.4)
	(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	Female	432 (47.9)	13 (48.1)	445 (47.9)	1.000
	Male	470 (52.1)	14 (51.9)	484 (52.1)
	(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	No	715 (79.3)	17 (63.0)	732 (78.8)	0.035
	Yes	166 (18.4)	10 (37.0)	176 (18.9)
	(Missing)	21 (2.3)	0 (0.0)	21 (2.3)

1.21 Dependent variable with any number of factor levels supported

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "extent.factor"
colon_s %>%
  dplyr::mutate(
    perfor.factor = ff_label(perfor.factor, "Perforated cancer")
  ) %>% 
  summary_factorlist(dependent, explanatory, p = TRUE, cont = "median", na_include = TRUE,
    column = TRUE, total_col = TRUE, add_dependent_label = TRUE, dependent_label_prefix = "") -> t

Extent of spread		Submucosa	Muscle	Serosa	Adjacent structures	Total	p
Age (years)	Median (IQR)	56.0 (50.0 to 64.0)	61.5 (56.0 to 70.0)	61.0 (53.0 to 69.0)	61.0 (53.5 to 66.0)	61.0 (53.0 to 69.0)	0.334
Age	<40 years	2 (9.5)	8 (7.5)	56 (7.4)	4 (9.3)	70 (7.5)	0.338
	40-59 years	12 (57.1)	32 (30.2)	285 (37.5)	15 (34.9)	344 (37.0)
	60+ years	7 (33.3)	66 (62.3)	418 (55.1)	24 (55.8)	515 (55.4)
	(Missing)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	Female	13 (61.9)	47 (44.3)	366 (48.2)	19 (44.2)	445 (47.9)	0.483
	Male	8 (38.1)	59 (55.7)	393 (51.8)	24 (55.8)	484 (52.1)
	(Missing)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	No	20 (95.2)	88 (83.0)	588 (77.5)	36 (83.7)	732 (78.8)	0.040
	Yes	1 (4.8)	13 (12.3)	157 (20.7)	5 (11.6)	176 (18.9)
	(Missing)	0 (0.0)	5 (4.7)	14 (1.8)	2 (4.7)	21 (2.3)

1.22 Missing data in the dependent

If you are careful to count totals and you know your data, you should realise when there is data missing from the dependent, e.g.:

explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "mort_5yr"
colon_s %>% 
    ff_glimpse(dependent, explanatory)
#> $Continuous
#>           label var_type   n missing_n missing_percent mean   sd  min
#> age Age (years)    <dbl> 929         0             0.0 59.8 11.9 18.0
#>     quartile_25 median quartile_75  max
#> age        53.0   61.0        69.0 85.0
#> 
#> $Categorical
#>                            label var_type   n missing_n missing_percent
#> mort_5yr        Mortality 5 year    <fct> 915        14             1.5
#> age.factor                   Age    <fct> 929         0             0.0
#> sex.factor                   Sex    <fct> 929         0             0.0
#> obstruct.factor      Obstruction    <fct> 908        21             2.3
#>                 levels_n                                               levels
#> mort_5yr               2                         "Alive", "Died", "(Missing)"
#> age.factor             3 "<40 years", "40-59 years", "60+ years", "(Missing)"
#> sex.factor             2                        "Female", "Male", "(Missing)"
#> obstruct.factor        2                             "No", "Yes", "(Missing)"
#>                 levels_count   levels_percent
#> mort_5yr        511, 404, 14 55.0, 43.5,  1.5
#> age.factor      70, 344, 515  7.5, 37.0, 55.4
#> sex.factor          445, 484           48, 52
#> obstruct.factor 732, 176, 21 78.8, 18.9,  2.3

To make sure, a warning is generated when data are dropped from the dependent:

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "mort_5yr"
colon_s %>%
  summary_factorlist(dependent, explanatory, p = TRUE, na_include = TRUE,
                                     total_col = TRUE,
                                     add_col_totals = TRUE, add_row_totals = TRUE) -> t
#> Note: dependent includes missing data. These are dropped.

label	Total N	Missing N	levels	Alive	Died	Total	p
Total N (%)				511 (55.8)	404 (44.2)	915
Age (years)	915 (100.0)	0	Mean (SD)	59.8 (11.4)	59.9 (12.5)	59.8 (11.9)	0.986
Age	915 (100.0)	0	<40 years	31 (6.1)	36 (8.9)	67 (7.3)	0.020
			40-59 years	208 (40.7)	131 (32.4)	339 (37.0)
			60+ years	272 (53.2)	237 (58.7)	509 (55.6)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	915 (100.0)	0	Female	243 (47.6)	194 (48.0)	437 (47.8)	0.941
			Male	268 (52.4)	210 (52.0)	478 (52.2)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	894 (97.7)	21	No	408 (79.8)	312 (77.2)	720 (78.7)	0.219
			Yes	89 (17.4)	85 (21.0)	174 (19.0)
			(Missing)	14 (2.7)	7 (1.7)	21 (2.3)

You may consider making the missing data explicit.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "mort_5yr"
colon_s %>%
    mutate(
        mort_5yr = forcats::fct_na_value_to_level(mort_5yr, level = "(Missing)")
    ) %>% 
  summary_factorlist(dependent, explanatory, p = TRUE, na_include = TRUE,
                                     total_col = TRUE,
                                     add_col_totals = TRUE, add_row_totals = TRUE) -> t

label	Total N	Missing N	levels	Alive	Died	(Missing)	Total	p
Total N (%)				511 (55.0)	404 (43.5)	14 (1.5)	929
Age (years)	929 (100.0)	0	Mean (SD)	59.8 (11.4)	59.9 (12.5)	53.9 (12.7)	59.8 (11.9)	0.185
Age	929 (100.0)	0	<40 years	31 (6.1)	36 (8.9)	3 (21.4)	70 (7.5)	0.018
			40-59 years	208 (40.7)	131 (32.4)	5 (35.7)	344 (37.0)
			60+ years	272 (53.2)	237 (58.7)	6 (42.9)	515 (55.4)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	929 (100.0)	0	Female	243 (47.6)	194 (48.0)	8 (57.1)	445 (47.9)	0.776
			Male	268 (52.4)	210 (52.0)	6 (42.9)	484 (52.1)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	908 (97.7)	21	No	408 (79.8)	312 (77.2)	12 (85.7)	732 (78.8)	0.373
			Yes	89 (17.4)	85 (21.0)	2 (14.3)	176 (18.9)
			(Missing)	14 (2.7)	7 (1.7)	0 (0.0)	21 (2.3)

1.23 Directly include missing data in dependent

Rather than making the data explicit in the dataset, you can use na_include_dependent = TRUE to do the same in summary_factorlist().

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "mort_5yr"
colon_s %>%
  summary_factorlist(dependent, explanatory, p = TRUE, 
                                     na_include = TRUE, na_include_dependent = TRUE,
                                     total_col = TRUE,
                                     add_col_totals = TRUE, add_row_totals = TRUE) -> t

label	Total N	Missing N	levels	Alive	Died	(Missing)	Total	p
Total N (%)				511 (55.0)	404 (43.5)	14 (1.5)	929
Age (years)	915 (100.0)	0	Mean (SD)	59.8 (11.4)	59.9 (12.5)	53.9 (12.7)	59.8 (11.9)	0.986
Age	915 (100.0)	0	<40 years	31 (6.1)	36 (8.9)	3 (21.4)	70 (7.5)	0.020
			40-59 years	208 (40.7)	131 (32.4)	5 (35.7)	344 (37.0)
			60+ years	272 (53.2)	237 (58.7)	6 (42.9)	515 (55.4)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	915 (100.0)	0	Female	243 (47.6)	194 (48.0)	8 (57.1)	445 (47.9)	0.941
			Male	268 (52.4)	210 (52.0)	6 (42.9)	484 (52.1)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	894 (97.7)	21	No	408 (79.8)	312 (77.2)	12 (85.7)	732 (78.8)	0.219
			Yes	89 (17.4)	85 (21.0)	2 (14.3)	176 (18.9)
			(Missing)	14 (2.7)	7 (1.7)	0 (0.0)	21 (2.3)

1.24 Summarise complete cases

You may wish to see summaries for complete cases across included variables. Rather than selecting including variables and drop_na(), you can pass na_complete_cases = TRUE to summary_factorlist() to do the same.

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "mort_5yr"
colon_s %>%
  summary_factorlist(dependent, explanatory, p = TRUE, 
                                     na_complete_cases = TRUE,
                                     total_col = TRUE,
                                     add_col_totals = TRUE, add_row_totals = TRUE) -> t
#> Note: dependent includes missing data. These are dropped.

label	Total N	Missing N	levels	Alive	Died	Total	p
Total N (%)				467 (55.7)	372 (44.3)	839
Age (years)	839 (91.7)	76	Mean (SD)	59.6 (11.5)	60.0 (12.3)	59.8 (11.9)	0.986
Age	839 (91.7)	76	<40 years	30 (6.4)	30 (8.1)	60 (7.2)	0.020
			40-59 years	190 (40.7)	122 (32.8)	312 (37.2)
			60+ years	247 (52.9)	220 (59.1)	467 (55.7)
Sex	839 (91.7)	76	Female	223 (47.8)	183 (49.2)	406 (48.4)	0.941
			Male	244 (52.2)	189 (50.8)	433 (51.6)
Obstruction	839 (91.7)	76	No	381 (81.6)	296 (79.6)	677 (80.7)	0.219
			Yes	86 (18.4)	76 (20.4)	162 (19.3)

1.25 Actively dropping missing data (and tidyverse functions that strip attributes)

You may wish to actively remove any rows with missing data, so you are explicit around which data are being used in models. Unfortunately some tidyverse functions silently remove variable attributes (labels). This is complained about then put right. But here is a workaround if it is happening with a variable you wish to use, such as tidyr::drop_na().

library(finalfit)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "mort_5yr"
vlabels = colon_s %>% 
    extract_variable_label()

colon_s %>%
    select(dependent, explanatory) %>% 
    tidyr::drop_na() %>% # Silently removes attributes
    ff_relabel(vlabels) %>% # Relabel
  summary_factorlist(dependent, explanatory, p = TRUE, na_include = TRUE,
                                     total_col = TRUE,
                                     add_col_totals = TRUE, add_row_totals = TRUE) -> t

label	Total N	Missing N	levels	Alive	Died	Total	p
Total N (%)				497 (55.6)	397 (44.4)	894
Age (years)	894 (100.0)	0	Mean (SD)	59.7 (11.5)	59.9 (12.5)	59.7 (11.9)	0.791
Age	894 (100.0)	0	<40 years	31 (6.2)	35 (8.8)	66 (7.4)	0.024
			40-59 years	203 (40.8)	129 (32.5)	332 (37.1)
			60+ years	263 (52.9)	233 (58.7)	496 (55.5)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Sex	894 (100.0)	0	Female	237 (47.7)	192 (48.4)	429 (48.0)	0.894
			Male	260 (52.3)	205 (51.6)	465 (52.0)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)
Obstruction	894 (100.0)	0	No	408 (82.1)	312 (78.6)	720 (80.5)	0.219
			Yes	89 (17.9)	85 (21.4)	174 (19.5)
			(Missing)	0 (0.0)	0 (0.0)	0 (0.0)

1.26 Explanatory variable defaults to factor when ≤5 distinct values

library(finalfit)

# Here, `extent` is a continuous variable with 4 distinct values. 
# Any continuous variable with 5 or fewer unique values is converted silently to factor 
# e.g.
explanatory = c("extent")
dependent = "mort_5yr"
colon_s %>%
  summary_factorlist(dependent, explanatory) -> t

label	levels	Alive	Died
extent	1	16 (3.1)	4 (1.0)
	2	78 (15.3)	25 (6.2)
	3	401 (78.5)	349 (86.4)
	4	16 (3.1)	26 (6.4)

1.27 Keep as continous variable when ≤5 distinct values

library(finalfit)

# Here, `extent` is a continuous variable with 4 distinct values. 
# Any continuous variable with 5 or fewer unique values is converted silently to factor 
# e.g.
explanatory = c("extent")
dependent = "mort_5yr"
colon_s %>%
  summary_factorlist(dependent, explanatory, 
                                     cont_cut = 0) -> t

label	levels	Alive	Died
extent	Mean (SD)	2.8 (0.5)	3.0 (0.4)

1.28 Stratified crosstables

I’ve been meaning to include support for table stratification for a while. I have delayed for a good reason. Perhaps the most straightforward way to implement stratificiation is with dplyr::group_by(). However, the non-standard evaluation required for multiple strata may confuse as it is not implemented else where in the package.

This translates to whether variable names are passed in quotes or not.

Here is a solution, which while not that pretty, is effective.

Note that tidyverse functions every so often start stripping labels/attributes. Hence the addition of the help function.

library(dplyr)
explanatory = c("age.factor", "sex.factor")
dependent = "perfor.factor"

# Pick option below
split = "rx.factor"
split = c("rx.factor", "node4.factor")

# Piped function to generate stratified crosstabs table
colon_s %>%
  group_by(!!! syms(split)) %>% # Looks awkward, but avoids unquoted var names
  group_modify(~ summary_factorlist(.x, dependent, explanatory)) %>%
  ff_stratify_helper(colon_s) -> t

Treatment	>4 positive nodes	levels	No	Yes
Obs	No	<40 years	14 (6.3)	0 (0.0)
Obs	No	40-59 years	89 (40.3)	3 (42.9)
Obs	No	60+ years	118 (53.4)	4 (57.1)
Obs	No	Female	101 (45.7)	3 (42.9)
Obs	No	Male	120 (54.3)	4 (57.1)
Obs	Yes	<40 years	10 (11.8)	1 (50.0)
Obs	Yes	40-59 years	31 (36.5)	1 (50.0)
Obs	Yes	60+ years	44 (51.8)	0 (0.0)
Obs	Yes	Female	44 (51.8)	1 (50.0)
Obs	Yes	Male	41 (48.2)	1 (50.0)
Lev	No	<40 years	14 (6.5)	0 (0.0)
Lev	No	40-59 years	78 (36.3)	3 (50.0)
Lev	No	60+ years	123 (57.2)	3 (50.0)
Lev	No	Female	89 (41.4)	2 (33.3)
Lev	No	Male	126 (58.6)	4 (66.7)
Lev	Yes	<40 years	4 (4.7)	1 (25.0)
Lev	Yes	40-59 years	33 (38.8)	1 (25.0)
Lev	Yes	60+ years	48 (56.5)	2 (50.0)
Lev	Yes	Female	39 (45.9)	3 (75.0)
Lev	Yes	Male	46 (54.1)	1 (25.0)
Lev+5FU	No	<40 years	15 (6.9)	0 (0.0)
Lev+5FU	No	40-59 years	72 (33.2)	2 (25.0)
Lev+5FU	No	60+ years	130 (59.9)	6 (75.0)
Lev+5FU	No	Female	115 (53.0)	4 (50.0)
Lev+5FU	No	Male	102 (47.0)	4 (50.0)
Lev+5FU	Yes	<40 years	11 (13.9)	0 (NaN)
Lev+5FU	Yes	40-59 years	31 (39.2)	0 (NaN)
Lev+5FU	Yes	60+ years	37 (46.8)	0 (NaN)
Lev+5FU	Yes	Female	44 (55.7)	0 (NaN)
Lev+5FU	Yes	Male	35 (44.3)	0 (NaN)

1.29 Digits / decimal places

explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
  summary_factorlist(dependent, explanatory, p = TRUE, digits = c(1,2,3,4,0)) -> t

label	levels	No	Yes	p
Age (years)	Mean (SD)	59.8 (11.91)	58.4 (13.30)	0.542
Age	<40 years	68 (7.5388)	2 (7.4074)	1.000
	40-59 years	334 (37.0288)	10 (37.0370)
	60+ years	500 (55.4324)	15 (55.5556)
Sex	Female	432 (47.8936)	13 (48.1481)	1.000
	Male	470 (52.1064)	14 (51.8519)
Obstruction	No	715 (81.1578)	17 (62.9630)	0.035
	Yes	166 (18.8422)	10 (37.0370)

1.30 Weighted tables

A simple weighting can be applied to tables. Explanatory continuous variables are multiplied by weights. Explanatory categorical variables are counted with a frequency weight (sum(weights)). This could be used with, say, inverse probability of treatment weightings (IPTW). The example uses random weights for demonstration purposes only. Hypothesis tests are not run on weighted data, p is set to FALSE.

explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
    mutate(my_weights = runif(929, 0, 1)) %>% # Random just to demonstrate
    summary_factorlist(dependent, explanatory, weights = "my_weights", digits = c(1, 1, 3, 1, 1))-> t

label	levels	No	Yes
Age (years)	Mean (SD)	59.5 (12.1)	59.3 (14.0)
Age	<40 years	37.3 (8.2)	1.0 (7.9)
	40-59 years	174.4 (38.3)	4.1 (31.6)
	60+ years	243.9 (53.5)	7.7 (60.4)
Sex	Female	218.5 (47.9)	4.6 (35.7)
	Male	237.2 (52.1)	8.2 (64.3)
Obstruction	No	362.7 (81.9)	7.7 (59.8)
	Yes	80.2 (18.1)	5.1 (40.2)

2 Model tables with `finalfit()`

2.01 Default

Logistic regression first.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.57 (0.34-0.98, p=0.041)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	0.81 (0.48-1.36, p=0.426)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.98 (0.75-1.28, p=0.902)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.25 (0.90-1.76, p=0.186)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.12 (0.51-2.44, p=0.770)

2.02 Hide reference levels

Most appropriate when all explanatory variables are continuous or well-known binary variables, such as sex.

library(finalfit)
explanatory = c("age", "sex.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory, add_dependent_label = FALSE) %>% 
    ff_remove_ref() %>% 
    dependent_label(colon_s, dependent)-> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age (years)	Mean (SD)	59.8 (11.4)	59.9 (12.5)	1.00 (0.99-1.01, p=0.986)	1.00 (0.99-1.01, p=0.983)
Sex	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.98 (0.76-1.27, p=0.888)

2.03 Model metrics

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory, metrics = TRUE) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.57 (0.34-0.98, p=0.041)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	0.81 (0.48-1.36, p=0.426)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.98 (0.75-1.28, p=0.902)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.25 (0.90-1.76, p=0.186)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.12 (0.51-2.44, p=0.770)

Number in dataframe = 929, Number in model = 894, Missing = 35, AIC = 1230.7, C-statistic = 0.56, H&L = Chi-sq(8) 5.69 (p=0.682)

2.04 Model metrics can be applied to all supported base models

library(finalfit)
glm(mort_5yr ~ age.factor + sex.factor + obstruct.factor + perfor.factor, data = colon_s, family = "binomial") %>% 
    ff_metrics() -> t

Number in dataframe = 929, Number in model = 894, Missing = 35, AIC = 1230.7, C-statistic = 0.56, H&L = Chi-sq(8) 5.69 (p=0.682)

2.05 Reduced model

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
explanatory_multi = c("age.factor", "obstruct.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory, explanatory_multi) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.57 (0.34-0.98, p=0.041)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	0.81 (0.48-1.36, p=0.424)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	-
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.26 (0.90-1.76, p=0.176)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	-

2.06 Include all models

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
explanatory_multi = c("age.factor", "obstruct.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory, explanatory_multi, metrics = TRUE, keep_models = TRUE) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable full)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.57 (0.34-0.98, p=0.041)	0.57 (0.34-0.98, p=0.041)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	0.81 (0.48-1.36, p=0.426)	0.81 (0.48-1.36, p=0.424)
Sex	Female	243 (55.6)	194 (44.4)	-	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.98 (0.75-1.28, p=0.902)	-
Obstruction	No	408 (56.7)	312 (43.3)	-	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.25 (0.90-1.76, p=0.186)	1.26 (0.90-1.76, p=0.176)
Perforation	No	497 (56.0)	391 (44.0)	-	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.12 (0.51-2.44, p=0.770)	-

Number in dataframe = 929, Number in model = 894, Missing = 35, AIC = 1230.7, C-statistic = 0.56, H&L = Chi-sq(8) 5.69 (p=0.682)

Number in dataframe = 929, Number in model = 894, Missing = 35, AIC = 1226.8, C-statistic = 0.555, H&L = Chi-sq(8) 0.06 (p=1.000)

2.06 Interactions

Interactions can be specified in the normal way. Formatting the output is trickier. At the moment, we have left the default model output. This can be adjusted as necessary.

library(finalfit)
explanatory = c("age.factor*sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.65 (0.32-1.34, p=0.241)	0.66 (0.32-1.36, p=0.258)
	60+ years	272 (53.4)	237 (46.6)	0.80 (0.40-1.61, p=0.529)	0.85 (0.42-1.71, p=0.647)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	1.24 (0.47-3.30, p=0.665)	1.17 (0.44-3.15, p=0.752)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.26 (0.90-1.76, p=0.182)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.11 (0.50-2.41, p=0.795)
age.factor40-59 years:sex.factorMale	Interaction	-	-	0.68 (0.23-1.97, p=0.479)	0.74 (0.25-2.18, p=0.588)
age.factor60+ years:sex.factorMale	Interaction	-	-	0.86 (0.30-2.39, p=0.766)	0.89 (0.31-2.51, p=0.822)

2.07 Interactions: create interaction variable with two factors

library(finalfit)
#explanatory = c("age.factor*sex.factor", "obstruct.factor", "perfor.factor")
explanatory = c("obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    ff_interaction(age.factor, sex.factor) %>% 
    finalfit(dependent, c(explanatory, "age.factor_sex.factor")) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.26 (0.90-1.76, p=0.182)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.11 (0.50-2.41, p=0.795)
Age:Sex	<40 years_Female	18 (48.6)	19 (51.4)	-	-
	<40 years_Male	13 (43.3)	17 (56.7)	1.24 (0.47-3.30, p=0.665)	1.17 (0.44-3.15, p=0.752)
	40-59 years_Female	96 (59.3)	66 (40.7)	0.65 (0.32-1.34, p=0.241)	0.66 (0.32-1.36, p=0.258)
	40-59 years_Male	112 (63.3)	65 (36.7)	0.55 (0.27-1.12, p=0.100)	0.57 (0.28-1.18, p=0.129)
	60+ years_Female	129 (54.2)	109 (45.8)	0.80 (0.40-1.61, p=0.529)	0.85 (0.42-1.71, p=0.647)
	60+ years_Male	143 (52.8)	128 (47.2)	0.85 (0.42-1.69, p=0.638)	0.88 (0.44-1.77, p=0.725)

2.08 Dependent name

The dependent name cannot be specified directly intentionally. This is to prevent errors when copying code. Re-label using ff_label(). The dependent prefix and suffix can also be altered.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    dplyr::mutate(
        mort_5yr = ff_label(mort_5yr, "5-year mortality")
    ) %>% 
    finalfit(dependent, explanatory, dependent_label_prefix = "",
                     dependent_label_suffix = " (full model)") -> t

5-year mortality (full model)		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.57 (0.34-0.98, p=0.041)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	0.81 (0.48-1.36, p=0.426)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.98 (0.75-1.28, p=0.902)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.25 (0.90-1.76, p=0.186)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.12 (0.51-2.44, p=0.770)

2.09 Estimate name

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory, estimate_name = "Odds ratio") -> t

Dependent: Mortality 5 year		Alive	Died	Odds ratio (univariable)	Odds ratio (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.57 (0.34-0.98, p=0.041)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	0.81 (0.48-1.36, p=0.426)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.98 (0.75-1.28, p=0.902)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.25 (0.90-1.76, p=0.186)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.12 (0.51-2.44, p=0.770)

2.10 Digits / decimal places

Number of digits to round to regression results. (1) estimate, (2) confidence interval limits, (3) p-value. Default is c(2,2,3). Trailing zeros are preserved. Number of decimal places for counts and mean (sd) / median (IQR) not currently supported. Defaults are senisble :)

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory, digits = c(3,3,4)) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.542 (0.319-0.918, p=0.0230)	0.574 (0.335-0.978, p=0.0412)
	60+ years	272 (53.4)	237 (46.6)	0.750 (0.448-1.250, p=0.2704)	0.810 (0.481-1.360, p=0.4261)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.981 (0.756-1.275, p=0.8886)	0.983 (0.754-1.283, p=0.9023)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.249 (0.896-1.741, p=0.1892)	1.255 (0.896-1.757, p=0.1859)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.180 (0.542-2.553, p=0.6716)	1.122 (0.512-2.442, p=0.7699)

2.11 Confidence interval type

One of c("profile", "default") for GLM models (confint.glm()). Note, a little awkwardly, the ‘default’ setting is profile, rather than default. Profile levels are probably a little more accurate. Only go to default if taking a significant length of time for profile, i.e. data is greater than hundreds of thousands of lines.

For glmer/lmer models (confint.merMod()), c("profile", "Wald", "boot"). Not implemented for lm(), coxph() or coxphlist, which use default.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory, confint_type = "default") -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.57 (0.34-0.98, p=0.041)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	0.81 (0.48-1.36, p=0.426)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.98 (0.75-1.28, p=0.902)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.25 (0.90-1.76, p=0.186)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.55-2.54, p=0.672)	1.12 (0.52-2.43, p=0.770)

2.12 Confidence interval level

Probably never change this :) Note, the p-value is intentionally not included for confidence levels other than 95% to avoid confusion.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory, confint_level = 0.90) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.35-0.84)	0.57 (0.37-0.90)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.49-1.15)	0.81 (0.52-1.25)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.79-1.22)	0.98 (0.79-1.23)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.95-1.65)	1.25 (0.95-1.66)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.62-2.25)	1.12 (0.58-2.15)

2.13 Confidence interval separation

Some like to avoid the hyphen so as not to confuse with minus sign. Obviously not an issue in logistic regression.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory, confint_sep = " to ") -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32 to 0.92, p=0.023)	0.57 (0.34 to 0.98, p=0.041)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45 to 1.25, p=0.270)	0.81 (0.48 to 1.36, p=0.426)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76 to 1.27, p=0.889)	0.98 (0.75 to 1.28, p=0.902)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90 to 1.74, p=0.189)	1.25 (0.90 to 1.76, p=0.186)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54 to 2.55, p=0.672)	1.12 (0.51 to 2.44, p=0.770)

2.14 Robust standard errors / confidence intervals

explanatory = c("age", "sex.factor")
dependent = 'mort_5yr'

# Standard finalfit regression table
t1 = colon_s %>%
  finalfit(dependent, explanatory, keep_fit_id = TRUE)

# GLM with Stata-like robust standard errors
t2 = colon_s %>% 
  glmmulti(dependent, explanatory) %>% 
  lmtest::coeftest(., vcov = sandwich::vcovHC(., "HC1")) %>% 
  broom::tidy(conf.int = TRUE) %>% 
  remove_intercept() %>% 
  select(term, estimate, conf.low, conf.high, p.value) %>% 
  mutate(across(c(estimate, conf.low, conf.high), exp)) %>% # or mutate_at(vars())
  as.data.frame() %>% 
  condense_fit(estimate_name = "OR (multivariable robust SE)")

ff_merge(t1, t2, last_merge = TRUE)

#> Note: dependent includes missing data. These are dropped.
#> Waiting for profiling to be done...
#> Waiting for profiling to be done...
#> Waiting for profiling to be done...

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)	OR (multivariable robust SE)
Age (years)	Mean (SD)	59.8 (11.4)	59.9 (12.5)	1.00 (0.99-1.01, p=0.986)	1.00 (0.99-1.01, p=0.983)	1.00 (0.99-1.01, p=0.984)
Sex	Female	243 (55.6)	194 (44.4)	-	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.98 (0.76-1.27, p=0.888)	0.98 (0.76-1.27, p=0.888)

2.15 Remove p-value

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
colon_s %>%
    finalfit(dependent, explanatory) %>% 
    ff_remove_p() -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92)	0.57 (0.34-0.98)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25)	0.81 (0.48-1.36)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27)	0.98 (0.75-1.28)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74)	1.25 (0.90-1.76)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55)	1.12 (0.51-2.44)

2.16 Mixed effects random-intercept model

At its simplest, a random-intercept model can be specified using a single quoted variable, e.g. random_effect = "hospital". This is equivalent to random_effect = "(1 | hospital)". Alternatively you can provide the full specification including parenthesis, e.g. random_effect = "(1 | hospital) + (1 | country)".

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
random_effect = "hospital"
colon_s %>%
    finalfit(dependent, explanatory, random_effect = random_effect,
                     dependent_label_suffix = " (random intercept)") -> t

Dependent: Mortality 5 year (random intercept)		Alive	Died	OR (univariable)	OR (multilevel)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.75 (0.39-1.44, p=0.382)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	1.03 (0.55-1.96, p=0.916)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.80 (0.58-1.11, p=0.180)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.23 (0.82-1.83, p=0.320)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.03 (0.43-2.51, p=0.940)

2.16b Mixed effects random-intercept model with univariable estimates including random effects

Some recently asked about this and it is a good question. Here is an approach.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
random_effect = "hospital"

colon_s %>%
    finalfit(dependent, explanatory, random_effect = random_effect, keep_fit_id = TRUE) %>% 
    ff_merge(
        explanatory %>% 
            purrr::map_df(~ glmmixed(colon_s, dependent, .x, random_effect = random_effect) %>% 
                                            fit2df(estimate_suffix = " (univariable with RE)")), 
        last_merge = TRUE
    ) %>% 
    dplyr::relocate(7, .before = 6) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (univariable with RE)	OR (multilevel)
Age	<40 years	31 (46.3)	36 (53.7)	-	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.68 (0.35-1.30, p=0.239)	0.75 (0.39-1.44, p=0.382)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	0.93 (0.49-1.74, p=0.809)	1.03 (0.55-1.96, p=0.916)
Sex	Female	243 (55.6)	194 (44.4)	-	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.81 (0.59-1.12, p=0.201)	0.80 (0.58-1.11, p=0.180)
Obstruction	No	408 (56.7)	312 (43.3)	-	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.23 (0.83-1.83, p=0.310)	1.23 (0.82-1.83, p=0.320)
Perforation	No	497 (56.0)	391 (44.0)	-	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.07 (0.44-2.57, p=0.888)	1.03 (0.43-2.51, p=0.940)

2.17 Mixed effects random-slope model

In the example below, allow the effect of age on outcome to vary by hospital. Note, this specification must have parentheses included.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
random_effect = "(age.factor | hospital)"
colon_s %>%
    finalfit(dependent, explanatory, random_effect = random_effect,
                     dependent_label_suffix = " (random slope: age)") -> t

Dependent: Mortality 5 year (random slope: age)		Alive	Died	OR (univariable)	OR (multilevel)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.54 (0.32-0.92, p=0.023)	0.81 (0.37-1.81, p=0.611)
	60+ years	272 (53.4)	237 (46.6)	0.75 (0.45-1.25, p=0.270)	1.08 (0.54-2.20, p=0.822)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	0.98 (0.76-1.27, p=0.889)	0.80 (0.58-1.11, p=0.179)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.24 (0.83-1.85, p=0.298)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.02 (0.42-2.48, p=0.967)

2.18 Mixed effects random-slope model directly from `lme4`

Clearly, as models get more complex, parameters such as random effect group variances may require to be extracted directly from model outputs.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
random_effect = "(age.factor | hospital)"
colon_s %>% 
    glmmixed(dependent, explanatory, random_effect = random_effect) %>% 
    broom.mixed::tidy() -> t

2.19 Exclude all missing data in final model from univariable analyses

This can be useful if you want the numbers in the final table to match the final multivariable model. However, be careful to include a full explanation of this in the methods and the reason for exluding the missing data.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = 'mort_5yr'
colon_s %>%
    dplyr::select(explanatory, dependent) %>%
    tidyr::drop_na() %>%
    finalfit(dependent, explanatory) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (47.0)	35 (53.0)	-	-
	40-59 years	203 (61.1)	129 (38.9)	0.56 (0.33-0.96, p=0.034)	0.57 (0.34-0.98, p=0.041)
	60+ years	263 (53.0)	233 (47.0)	0.78 (0.47-1.31, p=0.356)	0.81 (0.48-1.36, p=0.426)
Sex	Female	237 (55.2)	192 (44.8)	-	-
	Male	260 (55.9)	205 (44.1)	0.97 (0.75-1.27, p=0.841)	0.98 (0.75-1.28, p=0.902)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.25 (0.90-1.76, p=0.186)
Perforation	No	483 (55.7)	384 (44.3)	-	-
	Yes	14 (51.9)	13 (48.1)	1.17 (0.54-2.53, p=0.691)	1.12 (0.51-2.44, p=0.770)

2.20 Linear regression

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = 'nodes'
colon_s %>%
    finalfit(dependent, explanatory) -> t

Dependent: nodes		unit	value	Coefficient (univariable)	Coefficient (multivariable)
Age	<40 years	Mean (sd)	4.7 (4.5)	-	-
	40-59 years	Mean (sd)	3.6 (3.3)	-1.14 (-2.08 to -0.21, p=0.016)	-1.21 (-2.16 to -0.26, p=0.012)
	60+ years	Mean (sd)	3.6 (3.6)	-1.19 (-2.10 to -0.28, p=0.010)	-1.25 (-2.18 to -0.33, p=0.008)
Sex	Female	Mean (sd)	3.7 (3.6)	-	-
	Male	Mean (sd)	3.6 (3.6)	-0.14 (-0.60 to 0.33, p=0.565)	-0.07 (-0.54 to 0.40, p=0.779)
Obstruction	No	Mean (sd)	3.7 (3.7)	-	-
	Yes	Mean (sd)	3.5 (3.2)	-0.24 (-0.83 to 0.36, p=0.435)	-0.31 (-0.91 to 0.29, p=0.313)
Perforation	No	Mean (sd)	3.7 (3.6)	-	-
	Yes	Mean (sd)	3.9 (2.8)	0.24 (-1.13 to 1.61, p=0.735)	0.28 (-1.09 to 1.66, p=0.686)

2.21 Mixed effects random-intercept linear regression

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "nodes"
random_effect = "hospital"
colon_s %>%
    finalfit(dependent, explanatory, random_effect = random_effect,
                     dependent_label_suffix = " (random intercept)") -> t

Dependent: nodes (random intercept)		unit	value	Coefficient (univariable)	Coefficient (multilevel)
Age	<40 years	Mean (sd)	4.7 (4.5)	-	-
	40-59 years	Mean (sd)	3.6 (3.3)	-1.14 (-2.08 to -0.21, p=0.016)	-0.79 (-1.65 to 0.07, p=0.035)
	60+ years	Mean (sd)	3.6 (3.6)	-1.19 (-2.10 to -0.28, p=0.010)	-0.98 (-1.81 to -0.14, p=0.011)
Sex	Female	Mean (sd)	3.7 (3.6)	-	-
	Male	Mean (sd)	3.6 (3.6)	-0.14 (-0.60 to 0.33, p=0.565)	-0.19 (-0.62 to 0.24, p=0.195)
Obstruction	No	Mean (sd)	3.7 (3.7)	-	-
	Yes	Mean (sd)	3.5 (3.2)	-0.24 (-0.83 to 0.36, p=0.435)	-0.37 (-0.92 to 0.17, p=0.091)
Perforation	No	Mean (sd)	3.7 (3.6)	-	-
	Yes	Mean (sd)	3.9 (2.8)	0.24 (-1.13 to 1.61, p=0.735)	0.23 (-1.01 to 1.48, p=0.357)

2.22 Mixed effects random-slope linear regression

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "nodes"
random_effect = "(age.factor | hospital)"
colon_s %>%
    finalfit(dependent, explanatory, random_effect = random_effect,
                     dependent_label_suffix = " (random slope: age)") -> t

Dependent: nodes (random slope: age)		unit	value	Coefficient (univariable)	Coefficient (multilevel)
Age	<40 years	Mean (sd)	4.7 (4.5)	-	-
	40-59 years	Mean (sd)	3.6 (3.3)	-1.14 (-2.08 to -0.21, p=0.016)	-0.76 (-1.73 to 0.22, p=0.065)
	60+ years	Mean (sd)	3.6 (3.6)	-1.19 (-2.10 to -0.28, p=0.010)	-0.93 (-1.77 to -0.08, p=0.016)
Sex	Female	Mean (sd)	3.7 (3.6)	-	-
	Male	Mean (sd)	3.6 (3.6)	-0.14 (-0.60 to 0.33, p=0.565)	-0.19 (-0.62 to 0.24, p=0.196)
Obstruction	No	Mean (sd)	3.7 (3.7)	-	-
	Yes	Mean (sd)	3.5 (3.2)	-0.24 (-0.83 to 0.36, p=0.435)	-0.34 (-0.88 to 0.21, p=0.112)
Perforation	No	Mean (sd)	3.7 (3.6)	-	-
	Yes	Mean (sd)	3.9 (2.8)	0.24 (-1.13 to 1.61, p=0.735)	0.20 (-1.05 to 1.45, p=0.377)

2.23 Cox proportional hazards model (survival / time to event)

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "Surv(time, status)"
colon_s %>%
    finalfit(dependent, explanatory) -> t

Dependent: Surv(time, status)		all	HR (univariable)	HR (multivariable)
Age	<40 years	70 (7.5)	-	-
	40-59 years	344 (37.0)	0.76 (0.53-1.09, p=0.132)	0.79 (0.55-1.13, p=0.196)
	60+ years	515 (55.4)	0.93 (0.66-1.31, p=0.668)	0.98 (0.69-1.40, p=0.926)
Sex	Female	445 (47.9)	-	-
	Male	484 (52.1)	1.01 (0.84-1.22, p=0.888)	1.02 (0.85-1.23, p=0.812)
Obstruction	No	732 (80.6)	-	-
	Yes	176 (19.4)	1.29 (1.03-1.62, p=0.028)	1.30 (1.03-1.64, p=0.026)
Perforation	No	902 (97.1)	-	-
	Yes	27 (2.9)	1.17 (0.70-1.95, p=0.556)	1.08 (0.64-1.81, p=0.785)

2.24 Cox proportional hazards model: change dependent label

As above, the dependent label cannot be specfied directly in the model to avoid errors. However, in survival modelling the surivial object specification can be long or awkward. Therefore, here is the work around.

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "Surv(time, status)"
colon_s %>%
    finalfit(dependent, explanatory, add_dependent_label = FALSE) %>% 
    dplyr::rename("Overall survival" = label) %>% 
    dplyr::rename(" " = levels) -> t

Overall survival		all	HR (univariable)	HR (multivariable)
Age	<40 years	70 (7.5)	-	-
	40-59 years	344 (37.0)	0.76 (0.53-1.09, p=0.132)	0.79 (0.55-1.13, p=0.196)
	60+ years	515 (55.4)	0.93 (0.66-1.31, p=0.668)	0.98 (0.69-1.40, p=0.926)
Sex	Female	445 (47.9)	-	-
	Male	484 (52.1)	1.01 (0.84-1.22, p=0.888)	1.02 (0.85-1.23, p=0.812)
Obstruction	No	732 (80.6)	-	-
	Yes	176 (19.4)	1.29 (1.03-1.62, p=0.028)	1.30 (1.03-1.64, p=0.026)
Perforation	No	902 (97.1)	-	-
	Yes	27 (2.9)	1.17 (0.70-1.95, p=0.556)	1.08 (0.64-1.81, p=0.785)

3 Model tables manually using `ff_merge()`

3.1 Basic table

Note summary_factorlist() needs argument, fit_id = TRUE.

library(finalfit)
library(dplyr)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"

## Crosstable
colon_s %>%
    summary_factorlist(dependent, explanatory, fit_id=TRUE) -> table_1

## Univariable
colon_s %>%
    glmuni(dependent, explanatory) %>%
    fit2df(estimate_suffix=" (univariable)") -> table_2

## Merge

table_1 %>% 
    ff_merge(table_2) %>% 
    select(-c(fit_id, index)) %>% 
    dependent_label(colon_s, dependent)-> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)
Age	<40 years	31 (6.1)	36 (8.9)	-
	40-59 years	208 (40.7)	131 (32.4)	0.54 (0.32-0.92, p=0.023)
	60+ years	272 (53.2)	237 (58.7)	0.75 (0.45-1.25, p=0.270)
Sex	Female	243 (47.6)	194 (48.0)	-
	Male	268 (52.4)	210 (52.0)	0.98 (0.76-1.27, p=0.889)
Obstruction	No	408 (82.1)	312 (78.6)	-
	Yes	89 (17.9)	85 (21.4)	1.25 (0.90-1.74, p=0.189)
Perforation	No	497 (97.3)	391 (96.8)	-
	Yes	14 (2.7)	13 (3.2)	1.18 (0.54-2.55, p=0.672)

3.2 Complex table (all in single pipe)

library(finalfit)
library(dplyr)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
random_effect = "hospital"
dependent = "mort_5yr"

# All in one pipe

colon_s %>%
    ## Crosstable
    summary_factorlist(dependent, explanatory, fit_id=TRUE)  %>% 
    
    ## Add univariable
    ff_merge(
        glmuni(colon_s, dependent, explanatory) %>%
            fit2df(estimate_suffix=" (univariable)")
    ) %>% 
    
    ## Add multivariable
    ff_merge(
        glmmulti(colon_s, dependent, explanatory) %>%
            fit2df(estimate_suffix=" (multivariable)")
    ) %>% 
    
    ## Add mixed effects
    ff_merge(
        glmmixed(colon_s, dependent, explanatory, random_effect) %>%
            fit2df(estimate_suffix=" (multilevel)"),
        last_merge = TRUE
    ) %>% 
    dependent_label(colon_s, dependent) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)	OR (multilevel)
Age	<40 years	31 (6.1)	36 (8.9)	-	-	-
	40-59 years	208 (40.7)	131 (32.4)	0.54 (0.32-0.92, p=0.023)	0.57 (0.34-0.98, p=0.041)	0.75 (0.39-1.44, p=0.382)
	60+ years	272 (53.2)	237 (58.7)	0.75 (0.45-1.25, p=0.270)	0.81 (0.48-1.36, p=0.426)	1.03 (0.55-1.96, p=0.916)
Sex	Female	243 (47.6)	194 (48.0)	-	-	-
	Male	268 (52.4)	210 (52.0)	0.98 (0.76-1.27, p=0.889)	0.98 (0.75-1.28, p=0.902)	0.80 (0.58-1.11, p=0.180)
Obstruction	No	408 (82.1)	312 (78.6)	-	-	-
	Yes	89 (17.9)	85 (21.4)	1.25 (0.90-1.74, p=0.189)	1.25 (0.90-1.76, p=0.186)	1.23 (0.82-1.83, p=0.320)
Perforation	No	497 (97.3)	391 (96.8)	-	-	-
	Yes	14 (2.7)	13 (3.2)	1.18 (0.54-2.55, p=0.672)	1.12 (0.51-2.44, p=0.770)	1.03 (0.43-2.51, p=0.940)

3.3 Other GLM models

Poisson

library(finalfit)
library(dplyr)

## Dobson (1990) Page 93: Randomized Controlled Trial :
counts = c(18,17,15,20,10,20,25,13,12)
outcome = gl(3,1,9)
treatment = gl(3,3)
d.AD <- data.frame(treatment, outcome, counts)

dependent = "counts"
explanatory = c("outcome", "treatment")

fit_uni = d.AD %>% 
    glmuni(dependent, explanatory, family = poisson) %>% 
    fit2df(estimate_name = "Rate ratio (univariable)")

fit_multi = d.AD %>% 
    glmmulti(dependent, explanatory, family = poisson) %>% 
    fit2df(estimate_name = "Rate ratio (multivariable)")

# All in one pipe
d.AD %>%
    ## Crosstable
    summary_factorlist(dependent, explanatory, cont = "median", fit_id=TRUE)  %>% 
    
    ## Add univariable
    ff_merge(fit_uni, estimate_name = "Rate ratio") %>% 
    
    ## Add multivariable
    ff_merge(fit_multi, estimate_name = "Rate ratio",
                     last_merge = TRUE) %>% 
    dependent_label(d.AD, dependent) -> t

Dependent: counts		unit	value	Rate ratio (univariable)	Rate ratio (multivariable)
outcome	1	Median (IQR)	20.0 (19.0 to 22.5)	-	-
	2	Median (IQR)	13.0 (11.5 to 15.0)	0.63 (0.42-0.94, p=0.025)	0.63 (0.42-0.94, p=0.025)
	3	Median (IQR)	15.0 (13.5 to 17.5)	0.75 (0.51-1.09, p=0.128)	0.75 (0.51-1.09, p=0.128)
treatment	1	Median (IQR)	17.0 (16.0 to 17.5)	-	-
	2	Median (IQR)	20.0 (15.0 to 20.0)	1.00 (0.67-1.48, p=1.000)	1.00 (0.67-1.48, p=1.000)
	3	Median (IQR)	13.0 (12.5 to 19.0)	1.00 (0.67-1.48, p=1.000)	1.00 (0.67-1.48, p=1.000)

Gamma

library(finalfit)
library(dplyr)

# A Gamma example, from McCullagh & Nelder (1989, pp. 300-2)
clotting <- data.frame(
    u = c(5,10,15,20,30,40,60,80,100),
    lot1 = c(118,58,42,35,27,25,21,19,18),
    lot2 = c(69,35,26,21,18,16,13,12,12))

dependent = "lot1"
explanatory = "log(u)"

fit_uni = clotting %>% 
    glmuni(dependent, explanatory, family = Gamma) %>% 
    fit2df(estimate_name = "Coefficient", exp = FALSE, digits = c(3,3,4))

# All in one pipe
clotting %>%
    ## Crosstable
    summary_factorlist(dependent, explanatory, cont = "median", fit_id=TRUE)  %>% 
    
    ## Add fit
    ff_merge(fit_uni, last_merge = TRUE) %>% 
    dependent_label(colon_s, dependent) -> t

Dependent: lot1		unit	value	Coefficient
u	[5.0,100.0]	Median (IQR)	27.0 (21.0 to 42.0)	-
NA	NA	NA	NA	0.015 (0.015-0.016, p<0.0001)

3.4 Weighted regression

Note this is updated August 2022. Weights can specified from data and summary_factorlist() table is now weighted.

library(finalfit)
library(dplyr)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"

colon_s %>%
    mutate(myweights = runif(dim(colon_s)[1])) %>% # random just for example
    finalfit(dependent, explanatory, weights = "myweights") -> t

Dependent: Mortality 5 year (weighted)		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	13 (41.3)	19 (58.7)	-	-
	40-59 years	102 (62.6)	61 (37.4)	0.42 (0.19-0.90, p=0.027)	0.43 (0.20-0.93, p=0.033)
	60+ years	134 (54.5)	111 (45.5)	0.59 (0.27-1.23, p=0.162)	0.61 (0.28-1.29, p=0.201)
Sex	Female	125 (57.4)	93 (42.6)	-	-
	Male	124 (55.7)	98 (44.3)	1.07 (0.74-1.57, p=0.710)	1.08 (0.74-1.59, p=0.685)
Obstruction	No	196 (56.6)	151 (43.4)	-	-
	Yes	45 (53.7)	38 (46.3)	1.12 (0.69-1.81, p=0.639)	1.12 (0.69-1.83, p=0.641)
Perforation	No	242 (56.6)	185 (43.4)	-	-
	Yes	7 (54.1)	6 (45.9)	1.11 (0.36-3.32, p=0.852)	1.07 (0.34-3.22, p=0.909)

3.5 Using base R functions

Note ff_formula() convenience function to make multivariable formula (y ~ x1 + x2 + x3 etc.) from a dependent and explanatory vector of names.

library(finalfit)
library(dplyr)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"

# All in one pipe

colon_s %>%
    ## Crosstable
    summary_factorlist(dependent, explanatory, fit_id=TRUE)  %>% 
    
    ## Add univariable
    ff_merge(
        glmuni(colon_s, dependent, explanatory) %>%
            fit2df(estimate_suffix=" (univariable)")
    ) %>% 
    
    ## Add multivariable
    ff_merge(
        glm(
            ff_formula(dependent, explanatory), data = colon_s, family = "binomial", weights = NULL
        ) %>%
            fit2df(estimate_suffix=" (multivariable)"),
        last_merge = TRUE
    ) %>% 
    dependent_label(colon_s, dependent) -> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (6.1)	36 (8.9)	-	-
	40-59 years	208 (40.7)	131 (32.4)	0.54 (0.32-0.92, p=0.023)	0.57 (0.34-0.98, p=0.041)
	60+ years	272 (53.2)	237 (58.7)	0.75 (0.45-1.25, p=0.270)	0.81 (0.48-1.36, p=0.426)
Sex	Female	243 (47.6)	194 (48.0)	-	-
	Male	268 (52.4)	210 (52.0)	0.98 (0.76-1.27, p=0.889)	0.98 (0.75-1.28, p=0.902)
Obstruction	No	408 (82.1)	312 (78.6)	-	-
	Yes	89 (17.9)	85 (21.4)	1.25 (0.90-1.74, p=0.189)	1.25 (0.90-1.76, p=0.186)
Perforation	No	497 (97.3)	391 (96.8)	-	-
	Yes	14 (2.7)	13 (3.2)	1.18 (0.54-2.55, p=0.672)	1.12 (0.51-2.44, p=0.770)

3.6 Edit table rows

This can be done as any dataframe would be edited.

library(finalfit)
library(dplyr)
explanatory = c("age.factor*sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"

# Run model for term test
fit <- glm(
    ff_formula(dependent, explanatory), 
    data=colon_s, family = binomial
)

# Not run
#term_test <- survey::regTermTest(fit, "age.factor:sex.factor")

# Run final table with results of term test
colon_s %>%
    finalfit(dependent, explanatory) %>%
    rbind(c(
        "age.factor:sex.factor (overall)",
        "Interaction",
        "-",
        "-",
        "-",
        paste0("p = 0.775")
    ))-> t

Dependent: Mortality 5 year		Alive	Died	OR (univariable)	OR (multivariable)
Age	<40 years	31 (46.3)	36 (53.7)	-	-
	40-59 years	208 (61.4)	131 (38.6)	0.65 (0.32-1.34, p=0.241)	0.66 (0.32-1.36, p=0.258)
	60+ years	272 (53.4)	237 (46.6)	0.80 (0.40-1.61, p=0.529)	0.85 (0.42-1.71, p=0.647)
Sex	Female	243 (55.6)	194 (44.4)	-	-
	Male	268 (56.1)	210 (43.9)	1.24 (0.47-3.30, p=0.665)	1.17 (0.44-3.15, p=0.752)
Obstruction	No	408 (56.7)	312 (43.3)	-	-
	Yes	89 (51.1)	85 (48.9)	1.25 (0.90-1.74, p=0.189)	1.26 (0.90-1.76, p=0.182)
Perforation	No	497 (56.0)	391 (44.0)	-	-
	Yes	14 (51.9)	13 (48.1)	1.18 (0.54-2.55, p=0.672)	1.11 (0.50-2.41, p=0.795)
age.factor40-59 years:sex.factorMale	Interaction	-	-	0.68 (0.23-1.97, p=0.479)	0.74 (0.25-2.18, p=0.588)
age.factor60+ years:sex.factorMale	Interaction	-	-	0.86 (0.30-2.39, p=0.766)	0.89 (0.31-2.51, p=0.822)
age.factor:sex.factor (overall)	Interaction	-	-	-	p = 0.775

3.7 Base model + individual explanatory variables

This was an email enquiry about how to build on a base model. The example request was in a survival context.

This has been updated August 2019. We have left the original example of building the table from scratch as a comparison.

ff_permute() allows combinations of variables to be built on a base model. See options on help page to,

include the base model and/or the full model,
present the permuted variables at the top or bottom of the table,
produce separate model tables, or the default which is a single table.

library(dplyr)
mydata = colon_s
explanatory_base = c("age.factor", "sex.factor")
explanatory_permute = c("obstruct.factor", "perfor.factor", "node4.factor")
dependent = "Surv(time, status)"

mydata %>% 
    ff_permute(dependent, explanatory_base, explanatory_permute) %>% 
    rename("Overall survival" = `Dependent: Surv(time, status)`, # optional tidying
                 `n (%)` = "all") -> t

Overall survival		n (%)	HR (univariable)	HR (multivariable) 1	HR (multivariable) 2	HR (multivariable) 3	HR (multivariable) 4	HR (multivariable) 5
Age	<40 years	70 (7.5)	-	-	-	-	-	-
	40-59 years	344 (37.0)	0.76 (0.53-1.09, p=0.132)	0.76 (0.53-1.08, p=0.129)	0.79 (0.55-1.13, p=0.198)	0.76 (0.53-1.08, p=0.127)	0.85 (0.59-1.22, p=0.379)	0.90 (0.63-1.30, p=0.590)
	60+ years	515 (55.4)	0.93 (0.66-1.31, p=0.668)	0.93 (0.66-1.31, p=0.660)	0.98 (0.69-1.40, p=0.931)	0.92 (0.65-1.31, p=0.656)	1.09 (0.77-1.55, p=0.615)	1.19 (0.83-1.69, p=0.346)
Sex	Female	445 (47.9)	-	-	-	-	-	-
	Male	484 (52.1)	1.01 (0.84-1.22, p=0.888)	1.02 (0.85-1.23, p=0.847)	1.02 (0.85-1.24, p=0.803)	1.02 (0.85-1.22, p=0.854)	1.04 (0.87-1.26, p=0.647)	1.05 (0.87-1.27, p=0.597)
Obstruction	No	732 (80.6)	-	-	-	-	-	-
	Yes	176 (19.4)	1.29 (1.03-1.62, p=0.028)	-	1.31 (1.04-1.64, p=0.022)	-	-	1.35 (1.07-1.70, p=0.011)
Perforation	No	902 (97.1)	-	-	-	-	-	-
	Yes	27 (2.9)	1.17 (0.70-1.95, p=0.556)	-	-	1.18 (0.70-1.97, p=0.535)	-	1.16 (0.69-1.94, p=0.581)
>4 positive nodes	No	674 (72.6)	-	-	-	-	-	-
	Yes	255 (27.4)	2.60 (2.15-3.14, p<0.001)	-	-	-	2.64 (2.18-3.19, p<0.001)	2.68 (2.21-3.26, p<0.001)

library(finalfit)
library(dplyr)

mydata = colon_s
base_explanatory = c("age.factor", "sex.factor")
explanatory = c("obstruct.factor", "perfor.factor", "node4.factor")
dependent = "Surv(time, status)"

mydata %>%
    # Counts
    summary_factorlist(dependent, c(base_explanatory,
                                                                    explanatory),
                                         column = TRUE,
                                         fit_id = TRUE) %>% 
    
    # Univariable
    ff_merge(
        coxphuni(mydata, dependent, c(base_explanatory, explanatory)) %>% 
            fit2df(estimate_suffix = " (Univariable)")
    ) %>% 
    
    # Base
    ff_merge(
        coxphmulti(mydata, dependent, base_explanatory) %>% 
            fit2df(estimate_suffix = " (Base model)")
    ) %>% 
    
    # Model 1
    ff_merge(
        coxphmulti(mydata, dependent, c(base_explanatory, explanatory[1])) %>% 
            fit2df(estimate_suffix = " (Model 1)")
    ) %>% 
    
    # Model 2
    ff_merge(
        coxphmulti(mydata, dependent, c(base_explanatory, explanatory[2])) %>% 
            fit2df(estimate_suffix = " (Model 2)")
    ) %>% 
    
    # Model 3
    ff_merge(
        coxphmulti(mydata, dependent, c(base_explanatory, explanatory[3])) %>% 
            fit2df(estimate_suffix = " (Model 3)")
    ) %>% 
    
    # Full
    ff_merge(
        coxphmulti(mydata, dependent, c(base_explanatory, explanatory)) %>% 
            fit2df(estimate_suffix = " (Full)"),
        last_merge = TRUE
    ) %>% 
    
    # Tidy-up
    rename("Overall survival" = label) %>% 
    rename(" " = levels) %>% 
    rename(`n (%)` = all) -> t

Overall survival		n (%)	HR (Univariable)	HR (Base model)	HR (Model 1)	HR (Model 2)	HR (Model 3)	HR (Full)
Age	<40 years	70 (7.5)	-	-	-	-	-	-
	40-59 years	344 (37.0)	0.76 (0.53-1.09, p=0.132)	0.76 (0.53-1.08, p=0.129)	0.79 (0.55-1.13, p=0.198)	0.76 (0.53-1.08, p=0.127)	0.85 (0.59-1.22, p=0.379)	0.90 (0.63-1.30, p=0.590)
	60+ years	515 (55.4)	0.93 (0.66-1.31, p=0.668)	0.93 (0.66-1.31, p=0.660)	0.98 (0.69-1.40, p=0.931)	0.92 (0.65-1.31, p=0.656)	1.09 (0.77-1.55, p=0.615)	1.19 (0.83-1.69, p=0.346)
Sex	Female	445 (47.9)	-	-	-	-	-	-
	Male	484 (52.1)	1.01 (0.84-1.22, p=0.888)	1.02 (0.85-1.23, p=0.847)	1.02 (0.85-1.24, p=0.803)	1.02 (0.85-1.22, p=0.854)	1.04 (0.87-1.26, p=0.647)	1.05 (0.87-1.27, p=0.597)
Obstruction	No	732 (80.6)	-	-	-	-	-	-
	Yes	176 (19.4)	1.29 (1.03-1.62, p=0.028)	-	1.31 (1.04-1.64, p=0.022)	-	-	1.35 (1.07-1.70, p=0.011)
Perforation	No	902 (97.1)	-	-	-	-	-	-
	Yes	27 (2.9)	1.17 (0.70-1.95, p=0.556)	-	-	1.18 (0.70-1.97, p=0.535)	-	1.16 (0.69-1.94, p=0.581)
>4 positive nodes	No	674 (72.6)	-	-	-	-	-	-
	Yes	255 (27.4)	2.60 (2.15-3.14, p<0.001)	-	-	-	2.64 (2.18-3.19, p<0.001)	2.68 (2.21-3.26, p<0.001)

4 Support for complex survey structures via `library(survey)`

4.1 Linear regression

Examples taken from survey::svyglm() help page.

library(survey)
library(dplyr)

data(api)
dependent = "api00"
explanatory = c("ell", "meals", "mobility")

# Label data frame
apistrat = apistrat %>%
  mutate(
  api00 = ff_label(api00, "API in 2000 (api00)"),
  ell = ff_label(ell, "English language learners (percent)(ell)"),
  meals = ff_label(meals, "Meals eligible (percent)(meals)"),
  mobility = ff_label(mobility, "First year at the school (percent)(mobility)"),
  sch.wide = ff_label(sch.wide, "School-wide target met (sch.wide)")
  )

# Linear example
dependent = "api00"
explanatory = c("ell", "meals", "mobility")

# Stratified design
dstrat = svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)

# Univariable fit
fit_uni = dstrat %>%
  svyglmuni(dependent, explanatory) %>%
  fit2df(estimate_suffix = " (univariable)")

# Multivariable fit
fit_multi = dstrat %>%
  svyglmmulti(dependent, explanatory) %>%
  fit2df(estimate_suffix = " (multivariable)")

# Pipe together
apistrat %>%
  summary_factorlist(dependent, explanatory, fit_id = TRUE) %>%
  ff_merge(fit_uni) %>%
  ff_merge(fit_multi, last_merge = TRUE) %>%
  dependent_label(apistrat, dependent) -> t

Dependent: API in 2000 (api00)		unit	value	Coefficient (univariable)	Coefficient (multivariable)
English language learners (percent)(ell)	[0.0,84.0]	Mean (sd)	652.8 (121.0)	-3.73 (-4.36–3.10, p<0.001)	-0.48 (-1.25-0.29, p=0.222)
Meals eligible (percent)(meals)	[0.0,100.0]	Mean (sd)	652.8 (121.0)	-3.38 (-3.71–3.05, p<0.001)	-3.14 (-3.70–2.58, p<0.001)
First year at the school (percent)(mobility)	[1.0,99.0]	Mean (sd)	652.8 (121.0)	-1.43 (-3.31-0.46, p=0.137)	0.23 (-0.55-1.00, p=0.567)

4.2 Binomial example

Note model family needs specified and exponentiation set to TRUE if desired.

library(survey)
library(dplyr)

data(api)
dependent = "sch.wide"
explanatory = c("ell", "meals", "mobility")

# Label data frame
apistrat = apistrat %>%
  mutate(
  api00 = ff_label(api00, "API in 2000 (api00)"),
  ell = ff_label(ell, "English language learners (percent)(ell)"),
  meals = ff_label(meals, "Meals eligible (percent)(meals)"),
  mobility = ff_label(mobility, "First year at the school (percent)(mobility)"),
  sch.wide = ff_label(sch.wide, "School-wide target met (sch.wide)")
  )
  
# Univariable fit
fit_uni = dstrat %>%
  svyglmuni(dependent, explanatory, family = "quasibinomial") %>%
  fit2df(exp = TRUE, estimate_name = "OR", estimate_suffix = " (univariable)")

# Multivariable fit
fit_multi = dstrat %>%
  svyglmmulti(dependent, explanatory, family = "quasibinomial") %>%
  fit2df(exp = TRUE, estimate_name = "OR", estimate_suffix = " (multivariable)")

# Pipe together
apistrat %>%
  summary_factorlist(dependent, explanatory, fit_id = TRUE) %>%
  ff_merge(fit_uni) %>%
  ff_merge(fit_multi, last_merge = TRUE) %>%
  dependent_label(apistrat, dependent) -> t

Dependent: School-wide target met (sch.wide)		No	Yes	OR (univariable)	OR (multivariable)
English language learners (percent)(ell)	Mean (SD)	22.5 (19.3)	20.5 (20.0)	1.00 (0.98-1.01, p=0.715)	1.00 (0.97-1.02, p=0.851)
Meals eligible (percent)(meals)	Mean (SD)	46.0 (29.1)	44.7 (29.0)	1.00 (0.99-1.01, p=0.968)	1.00 (0.98-1.02, p=0.732)
First year at the school (percent)(mobility)	Mean (SD)	13.9 (8.6)	17.2 (13.0)	1.06 (1.00-1.12, p=0.049)	1.06 (1.00-1.13, p=0.058)

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.

All tables examples

Ewen Harrison

1 Cross tables

1.01 Default

1.02 Add or edit variable labels

1.03 P-value for hypothesis test

1.04 With Fisher’s exact test

1.05 Parametric explanatory variables

1.06 Non-parametric explanatory variables

1.07 Select specific non-parametric variables

1.08 Missing values for the explanatory variables

1.09 Pass missing values to statistical tests

1.10 Row proportions (rather than column)

1.11 Total column

1.12 Row totals with missing

1.13 Row totals without missing

1.14 Row totals with user-specified column names

1.15 Order a variable by total

1.17 Add column totals

1.18 Add column totals without proportion.

1.19 Add column totals with user-specified row name and prefix.

1.20 Label with dependent name

1.21 Dependent variable with any number of factor levels supported

1.22 Missing data in the dependent

1.23 Directly include missing data in dependent

1.24 Summarise complete cases

1.25 Actively dropping missing data (and tidyverse functions that strip attributes)

1.26 Explanatory variable defaults to factor when ≤5 distinct values

1.27 Keep as continous variable when ≤5 distinct values

1.28 Stratified crosstables

1.29 Digits / decimal places

1.30 Weighted tables

2 Model tables with finalfit()

2.01 Default

2.02 Hide reference levels

2.03 Model metrics

2.04 Model metrics can be applied to all supported base models

2.05 Reduced model

2.06 Include all models

2.06 Interactions

2.07 Interactions: create interaction variable with two factors

2.08 Dependent name

2.09 Estimate name

2.10 Digits / decimal places

2.11 Confidence interval type

2.12 Confidence interval level

2.13 Confidence interval separation

2.14 Robust standard errors / confidence intervals

2.15 Remove p-value

2.16 Mixed effects random-intercept model

2.16b Mixed effects random-intercept model with univariable estimates including random effects

2.17 Mixed effects random-slope model

2.18 Mixed effects random-slope model directly from lme4

2.19 Exclude all missing data in final model from univariable analyses

2.20 Linear regression

2.21 Mixed effects random-intercept linear regression

2.22 Mixed effects random-slope linear regression

2.23 Cox proportional hazards model (survival / time to event)

2.24 Cox proportional hazards model: change dependent label

3 Model tables manually using ff_merge()

3.1 Basic table

3.2 Complex table (all in single pipe)

3.3 Other GLM models

Poisson

Gamma

3.4 Weighted regression

3.5 Using base R functions

3.6 Edit table rows

3.7 Base model + individual explanatory variables

4 Support for complex survey structures via library(survey)

4.1 Linear regression

4.2 Binomial example

1.20 Label with `dependent` name

2 Model tables with `finalfit()`

2.18 Mixed effects random-slope model directly from `lme4`

3 Model tables manually using `ff_merge()`

4 Support for complex survey structures via `library(survey)`