# All tables examples

install.packages("finalfit")

## 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 (maximum of 5 levels) - 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
#> Warning in chisq.test(tab, correct = FALSE): Chi-squared approximation may
#> be incorrect
label levels No Yes
Age (years) Mean (SD) 59.8 (11.9) 58.4 (13.3)
Age <40 years 68 (97.1) 2 (2.9)
40-59 years 334 (97.1) 10 (2.9)
60+ years 500 (97.1) 15 (2.9)
Sex Female 432 (97.1) 13 (2.9)
Male 470 (97.1) 14 (2.9)
Obstruction No 715 (97.7) 17 (2.3)
Yes 166 (94.3) 10 (5.7)

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)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#>     filter, lag
#> The following objects are masked from 'package:base':
#>
#>     intersect, setdiff, setequal, union
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 (97.1) 2 (2.9)
40-59 years 334 (97.1) 10 (2.9)
60+ years 500 (97.1) 15 (2.9)
Gender Female 432 (97.1) 13 (2.9)
Male 470 (97.1) 14 (2.9)
Obstruction No 715 (97.7) 17 (2.3)
Yes 166 (94.3) 10 (5.7)

### 1.03 P-value for hypothesis test

Chi-squared for categorical, Kruskal-Wallis/Mann-Whitney for continuous

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.578
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)
Sex Female 432 (97.1) 13 (2.9) 0.979
Male 470 (97.1) 14 (2.9)
Obstruction No 715 (97.7) 17 (2.3) 0.018
Yes 166 (94.3) 10 (5.7)

### 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, catTest = catTestfisher) -> t
label levels No Yes p
Age (years) Mean (SD) 59.8 (11.9) 58.4 (13.3) 0.578
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)
Sex Female 432 (97.1) 13 (2.9) 1.000
Male 470 (97.1) 14 (2.9)
Obstruction No 715 (97.7) 17 (2.3) 0.026
Yes 166 (94.3) 10 (5.7)

### 1.05 Median (interquartile range) instead of mean (standard deviation)

… for continuous variables.

library(finalfit)
explanatory = c("age", "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 (16.0) 60.0 (18.0) 0.578
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)
Sex Female 432 (97.1) 13 (2.9) 0.979
Male 470 (97.1) 14 (2.9)
Obstruction No 715 (97.7) 17 (2.3) 0.018
Yes 166 (94.3) 10 (5.7)

### 1.06 Missing values for the explanatory variables

Always do this 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, cont = "median", na_include = TRUE) -> t
label levels No Yes p
Age (years) Median (IQR) 61.0 (16.0) 60.0 (18.0) 0.578
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)
Sex Female 432 (97.1) 13 (2.9) 0.979
Male 470 (97.1) 14 (2.9)
Obstruction No 715 (97.7) 17 (2.3) 0.042
Yes 166 (94.3) 10 (5.7)
Missing 21 (100.0) 0 (0.0)

### 1.07 Column proportions (rather than row)

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) -> t
label levels No Yes p
Age (years) Median (IQR) 61.0 (16.0) 60.0 (18.0) 0.578
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) 0.979
Male 470 (52.1) 14 (51.9)
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.08 Total column

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 (16.0) 60.0 (18.0) 61.0 (16.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)
Sex Female 432 (47.9) 13 (48.1) 445 (47.9) 0.979
Male 470 (52.1) 14 (51.9) 484 (52.1)
Obstruction No 715 (79.3) 17 (63.0) 732 (78.8) 0.042
Yes 166 (18.4) 10 (37.0) 176 (18.9)
Missing 21 (2.3) 0 (0.0) 21 (2.3)

### 1.09 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)

### 1.10 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 (16.0) 60.0 (18.0) 61.0 (16.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)
Sex Female 432 (47.9) 13 (48.1) 445 (47.9) 0.979
Male 470 (52.1) 14 (51.9) 484 (52.1)
Obstruction No 715 (79.3) 17 (63.0) 732 (78.8) 0.042
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 (16.0) 60.0 (18.0) 61.0 (16.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)
Sex Female 432 (47.9) 13 (48.1) 445 (47.9) 0.979
Male 470 (52.1) 14 (51.9) 484 (52.1)
Obstruction No 715 (79.3) 17 (63.0) 732 (78.8) 0.042
Yes 166 (18.4) 10 (37.0) 176 (18.9)
Missing 21 (2.3) 0 (0.0) 21 (2.3)

### 1.11 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
Age (years) Median (IQR) 61.0 (12.5) 61.5 (14.0) 61.0 (16.0) 56.0 (14.0) 61.0 (16.0) 0.334
Age <40 years 4 (9.3) 8 (7.5) 56 (7.4) 2 (9.5) 70 (7.5) 0.338
40-59 years 15 (34.9) 32 (30.2) 285 (37.5) 12 (57.1) 344 (37.0)
60+ years 24 (55.8) 66 (62.3) 418 (55.1) 7 (33.3) 515 (55.4)
Sex Female 19 (44.2) 47 (44.3) 366 (48.2) 13 (61.9) 445 (47.9) 0.483
Male 24 (55.8) 59 (55.7) 393 (51.8) 8 (38.1) 484 (52.1)
Obstruction No 36 (83.7) 88 (83.0) 588 (77.5) 20 (95.2) 732 (78.8) 0.037
Yes 5 (11.6) 13 (12.3) 157 (20.7) 1 (4.8) 176 (18.9)
Missing 2 (4.7) 5 (4.7) 14 (1.8) 0 (0.0) 21 (2.3)

### 1.12 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 (80.0) 4 (20.0)
2 78 (75.7) 25 (24.3)
3 401 (53.5) 349 (46.5)
4 16 (38.1) 26 (61.9)

### 1.13 Keep as continous variable when ≤5 distinct values

library(finalfit)
explanatory = c("extent")
dependent = "mort_5yr"
colon_s %>%
summary_factorlist(dependent, explanatory, cont_cut = 3) -> t
label levels Alive Died
extent Mean (SD) 2.8 (0.5) 3.0 (0.4)

### 1.14 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 (doesn’t work with group_by_). This translates to whether variable names are passed in quotes or not. Finally,. dplyr::do() is planned for deprecation, but there is no good alternative at the moment. Anyway, here is a solution, which while not that pretty, is very effective.

library(dplyr)
# Piped function to generate stratified crosstabs table
explanatory = c("age.factor", "sex.factor")
dependent = "rx.factor"

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

colon_s %>%
group_by(!!! syms(split)) %>% #Looks awkward, but this keeps quoted var names (rather than unquoted)
do(
summary_factorlist(., dependent, explanatory, p = TRUE)
) %>%
data.frame() %>%
dependent_label(colon_s, dependent, prefix = "") %>%
colname2label(split) -> t
Perforation >4 positive nodes Treatment Obs Lev Lev.5FU p
No No Age <40 years 14 (32.6) 14 (32.6) 15 (34.9) 0.663
No No 40-59 years 89 (37.2) 78 (32.6) 72 (30.1)
No No 60+ years 118 (31.8) 123 (33.2) 130 (35.0)
No No Sex Female 101 (33.1) 89 (29.2) 115 (37.7) 0.050
No No Male 120 (34.5) 126 (36.2) 102 (29.3)
No Yes Age <40 years 10 (40.0) 4 (16.0) 11 (44.0) 0.322
No Yes 40-59 years 31 (32.6) 33 (34.7) 31 (32.6)
No Yes 60+ years 44 (34.1) 48 (37.2) 37 (28.7)
No Yes Sex Female 44 (34.6) 39 (30.7) 44 (34.6) 0.448
No Yes Male 41 (33.6) 46 (37.7) 35 (28.7)
Yes No Age <40 years 0 (NaN) 0 (NaN) 0 (NaN) 0.604
Yes No 40-59 years 3 (37.5) 3 (37.5) 2 (25.0)
Yes No 60+ years 4 (30.8) 3 (23.1) 6 (46.2)
Yes No Sex Female 3 (33.3) 2 (22.2) 4 (44.4) 0.823
Yes No Male 4 (33.3) 4 (33.3) 4 (33.3)
Yes Yes Age <40 years 1 (50.0) 1 (50.0) NA 0.472
Yes Yes 40-59 years 1 (50.0) 1 (50.0) NA
Yes Yes 60+ years 0 (0.0) 2 (100.0) NA
Yes Yes Sex Female 1 (25.0) 3 (75.0) NA 0.540
Yes Yes Male 1 (50.0) 1 (50.0) NA

## 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 (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)

### 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 (52.4) 210 (52.0) 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 (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)
 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 (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.424)
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) 1.26 (0.90-1.76, p=0.176)
Perforation No 497 (97.3) 391 (96.8) - -
Yes 14 (2.7) 13 (3.2) 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) OR (multivariable reduced)
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.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) 0.81 (0.48-1.36, p=0.424)
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) 1.26 (0.90-1.76, p=0.176)
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) -
 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 (6.1) 36 (8.9) - -
40-59 years 208 (40.7) 131 (32.4) 0.65 (0.32-1.34, p=0.241) 0.66 (0.32-1.36, p=0.258)
60+ years 272 (53.2) 237 (58.7) 0.80 (0.40-1.61, p=0.529) 0.85 (0.42-1.71, p=0.647)
Sex Female 243 (47.6) 194 (48.0) - -
Male 268 (52.4) 210 (52.0) 1.24 (0.47-3.30, p=0.665) 1.17 (0.44-3.15, p=0.752)
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.26 (0.90-1.76, p=0.182)
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.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 (82.1) 312 (78.6) - -
Yes 89 (17.9) 85 (21.4) 1.25 (0.90-1.74, p=0.189) 1.26 (0.90-1.76, p=0.182)
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.11 (0.50-2.41, p=0.795)
Age:Sex <40 years|Female 18 (3.5) 19 (4.7) - -
<40 years|Male 13 (2.5) 17 (4.2) 1.24 (0.47-3.30, p=0.665) 1.17 (0.44-3.15, p=0.752)
40-59 years|Female 96 (18.8) 66 (16.3) 0.65 (0.32-1.34, p=0.241) 0.66 (0.32-1.36, p=0.258)
40-59 years|Male 112 (21.9) 65 (16.1) 0.55 (0.27-1.12, p=0.100) 0.57 (0.28-1.18, p=0.129)
60+ years|Female 129 (25.2) 109 (27.0) 0.80 (0.40-1.61, p=0.529) 0.85 (0.42-1.71, p=0.647)
60+ years|Male 143 (28.0) 128 (31.7) 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 (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)

### 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 (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)

### 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 (6.1) 36 (8.9) - -
40-59 years 208 (40.7) 131 (32.4) 0.542 (0.319-0.918, p=0.0230) 0.574 (0.335-0.978, p=0.0412)
60+ years 272 (53.2) 237 (58.7) 0.750 (0.448-1.250, p=0.2704) 0.810 (0.481-1.360, p=0.4261)
Sex Female 243 (47.6) 194 (48.0) - -
Male 268 (52.4) 210 (52.0) 0.981 (0.756-1.275, p=0.8886) 0.983 (0.754-1.283, p=0.9023)
Obstruction No 408 (82.1) 312 (78.6) - -
Yes 89 (17.9) 85 (21.4) 1.249 (0.896-1.741, p=0.1892) 1.255 (0.896-1.757, p=0.1859)
Perforation No 497 (97.3) 391 (96.8) - -
Yes 14 (2.7) 13 (3.2) 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 (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.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 (6.1) 36 (8.9) - -
40-59 years 208 (40.7) 131 (32.4) 0.54 (0.35-0.84) 0.57 (0.37-0.90)
60+ years 272 (53.2) 237 (58.7) 0.75 (0.49-1.15) 0.81 (0.52-1.25)
Sex Female 243 (47.6) 194 (48.0) - -
Male 268 (52.4) 210 (52.0) 0.98 (0.79-1.22) 0.98 (0.79-1.23)
Obstruction No 408 (82.1) 312 (78.6) - -
Yes 89 (17.9) 85 (21.4) 1.25 (0.95-1.65) 1.25 (0.95-1.66)
Perforation No 497 (97.3) 391 (96.8) - -
Yes 14 (2.7) 13 (3.2) 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 (6.1) 36 (8.9) - -
40-59 years 208 (40.7) 131 (32.4) 0.54 (0.32 to 0.92, p=0.023) 0.57 (0.34 to 0.98, p=0.041)
60+ years 272 (53.2) 237 (58.7) 0.75 (0.45 to 1.25, p=0.270) 0.81 (0.48 to 1.36, p=0.426)
Sex Female 243 (47.6) 194 (48.0) - -
Male 268 (52.4) 210 (52.0) 0.98 (0.76 to 1.27, p=0.889) 0.98 (0.75 to 1.28, p=0.902)
Obstruction No 408 (82.1) 312 (78.6) - -
Yes 89 (17.9) 85 (21.4) 1.25 (0.90 to 1.74, p=0.189) 1.25 (0.90 to 1.76, p=0.186)
Perforation No 497 (97.3) 391 (96.8) - -
Yes 14 (2.7) 13 (3.2) 1.18 (0.54 to 2.55, p=0.672) 1.12 (0.51 to 2.44, p=0.770)

### 2.14 Mixed effects random-intercept model

At its simplest, a random-intercept model can be specified using a single quoted variable. In this example, it is the equivalent of quoting random_effect = "(1 | hospital)".

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 (6.1) 36 (8.9) - -
40-59 years 208 (40.7) 131 (32.4) 0.54 (0.32-0.92, p=0.023) 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) 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.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.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.03 (0.43-2.51, p=0.940)

### 2.15 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 (6.1) 36 (8.9) - -
40-59 years 208 (40.7) 131 (32.4) 0.54 (0.32-0.92, p=0.023) 0.81 (0.37-1.81, p=0.611)
60+ years 272 (53.2) 237 (58.7) 0.75 (0.45-1.25, p=0.270) 1.08 (0.54-2.20, p=0.822)
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.80 (0.58-1.11, p=0.179)
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.24 (0.83-1.85, p=0.298)
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.02 (0.42-2.48, p=0.967)

### 2.16 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 %>%
lme4::glmer(mort_5yr ~ age.factor + (age.factor | hospital), family = "binomial", data = .) %>%
broom::tidy() -> t
term estimate std.error statistic p.value group
(Intercept) -0.2537662 0.8983060 -0.2824941 0.7775646 fixed
age.factor40-59 years -0.3285638 0.3830047 -0.8578582 0.3909707 fixed
age.factor60+ years -0.0531730 0.3450263 -0.1541128 0.8775208 fixed
sd_(Intercept).hospital 1.8670680 NA NA NA hospital
sd_age.factor40-59 years.hospital 0.3382630 NA NA NA hospital
sd_age.factor60+ years.hospital 0.0826644 NA NA NA hospital
cor_(Intercept).age.factor40-59 years.hospital -0.9999999 NA NA NA hospital
cor_(Intercept).age.factor60+ years.hospital -0.9999997 NA NA NA hospital
cor_age.factor40-59 years.age.factor60+ years.hospital 0.9999998 NA NA NA hospital

### 2.17 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) %>%
na.omit() %>%
finalfit(dependent, explanatory) -> t
Dependent: Mortality 5 year Alive Died OR (univariable) OR (multivariable)
Age <40 years 31 (6.2) 35 (8.8) - -
40-59 years 203 (40.8) 129 (32.5) 0.56 (0.33-0.96, p=0.034) 0.57 (0.34-0.98, p=0.041)
60+ years 263 (52.9) 233 (58.7) 0.78 (0.47-1.31, p=0.356) 0.81 (0.48-1.36, p=0.426)
Sex Female 237 (47.7) 192 (48.4) - -
Male 260 (52.3) 205 (51.6) 0.97 (0.75-1.27, p=0.841) 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 483 (97.2) 384 (96.7) - -
Yes 14 (2.8) 13 (3.3) 1.17 (0.54-2.53, p=0.691) 1.12 (0.51-2.44, p=0.770)

### 2.18 Linear regression

library(finalfit)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = 'nodes'
colon_s %>%
finalfit(dependent, explanatory) -> t
Dependent: nodes Mean (sd) Coefficient (univariable) Coefficient (multivariable)
Age <40 years 4.7 (4.5) - -
40-59 years 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 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 3.7 (3.6) - -
Male 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 3.7 (3.7) - -
Yes 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 3.7 (3.6) - -
Yes 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.19 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) Mean (sd) Coefficient (univariable) Coefficient (multilevel)
Age <40 years 4.7 (4.5) - -
40-59 years 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 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 3.7 (3.6) - -
Male 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 3.7 (3.7) - -
Yes 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 3.7 (3.6) - -
Yes 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.20 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) Mean (sd) Coefficient (univariable) Coefficient (multilevel)
Age <40 years 4.7 (4.5) - -
40-59 years 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 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 3.7 (3.6) - -
Male 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 3.7 (3.7) - -
Yes 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 3.7 (3.6) - -
Yes 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.21 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.22 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 (46.3) 36 (53.7) -
40-59 years 208 (61.4) 131 (38.6) 0.54 (0.32-0.92, p=0.023)
60+ years 272 (53.4) 237 (46.6) 0.75 (0.45-1.25, p=0.270)
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)
Perforation No 497 (56.0) 391 (44.0) -
Yes 14 (51.9) 13 (48.1) 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)  %>%

ff_merge(
glmuni(colon_s, dependent, explanatory) %>%
fit2df(estimate_suffix=" (univariable)")
) %>%

ff_merge(
glmmulti(colon_s, dependent, explanatory) %>%
fit2df(estimate_suffix=" (multivariable)")
) %>%

ff_merge(
glmmixed(colon_s, dependent, explanatory, random_effect) %>%
fit2df(estimate_suffix=" (multilevel)")
) %>%
select(-c(fit_id, index)) %>%
dependent_label(colon_s, dependent) -> t
Dependent: Mortality 5 year Alive Died OR (univariable) OR (multivariable) 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.57 (0.34-0.98, p=0.041) 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.81 (0.48-1.36, p=0.426) 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.98 (0.75-1.28, p=0.902) 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.25 (0.90-1.76, p=0.186) 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.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)

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

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

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

# All in one pipe
## Crosstable
summary_factorlist(dependent, explanatory, cont = "median", fit_id=TRUE)  %>%

ff_merge(fit_uni, estimate_name = "Rate ratio") %>%

ff_merge(fit_multi, estimate_name = "Rate ratio") %>%

select(-c(fit_id, index)) %>%
dependent_label(d.AD, dependent) -> t
Dependent: counts Median (IQR) Rate ratio (univariable) Rate ratio (multivariable)
outcome 1 20.0 (19.0 to 22.5) - -
2 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 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 17.0 (16.0 to 17.5) - -
2 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 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)  %>%

ff_merge(fit_uni) %>%

select(-c(fit_id, index)) %>%
dependent_label(colon_s, dependent) -> t
Dependent: lot1 Median (IQR) Coefficient
log(u) [1.61,4.61] 27.0 (21.0 to 42.0) 0.015 (0.015-0.016, p<0.0001)

### 3.4 Weighted regression

library(finalfit)
library(dplyr)
explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = "mort_5yr"
weights = runif(dim(colon_s)[1]) # random just for example

# All in one pipe
colon_s %>%
## Crosstable
summary_factorlist(dependent, explanatory, fit_id=TRUE)  %>%

ff_merge(
glmuni(colon_s, dependent, explanatory, weights = weights, family = quasibinomial) %>%
fit2df(estimate_suffix=" (univariable)")
) %>%

ff_merge(
glmmulti(colon_s, dependent, explanatory, weights = weights, family = quasibinomial) %>%
fit2df(estimate_suffix=" (multivariable)")
) %>%
select(-c(fit_id, index)) %>%
dependent_label(colon_s, dependent) -> 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.51 (0.29-0.87, p=0.014) 0.54 (0.31-0.94, p=0.030)
60+ years 272 (53.4) 237 (46.6) 0.76 (0.45-1.28, p=0.305) 0.83 (0.49-1.41, p=0.493)
Sex Female 243 (55.6) 194 (44.4) - -
Male 268 (56.1) 210 (43.9) 0.91 (0.70-1.19, p=0.494) 0.92 (0.70-1.20, p=0.517)
Obstruction No 408 (56.7) 312 (43.3) - -
Yes 89 (51.1) 85 (48.9) 1.16 (0.82-1.62, p=0.401) 1.18 (0.83-1.66, p=0.354)
Perforation No 497 (56.0) 391 (44.0) - -
Yes 14 (51.9) 13 (48.1) 1.00 (0.46-2.13, p=0.993) 1.00 (0.46-2.14, p=0.994)

### 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)  %>%

ff_merge(
glmuni(colon_s, dependent, explanatory) %>%
fit2df(estimate_suffix=" (univariable)")
) %>%

ff_merge(
glm(
ff_formula(dependent, explanatory), data = colon_s, family = "binomial", weights = NULL
) %>%
fit2df(estimate_suffix=" (multivariable)")
) %>%

select(-c(fit_id, index)) %>%
dependent_label(colon_s, dependent) -> 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)

### 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 (6.1) 36 (8.9) - -
40-59 years 208 (40.7) 131 (32.4) 0.65 (0.32-1.34, p=0.241) 0.66 (0.32-1.36, p=0.258)
60+ years 272 (53.2) 237 (58.7) 0.80 (0.40-1.61, p=0.529) 0.85 (0.42-1.71, p=0.647)
Sex Female 243 (47.6) 194 (48.0) - -
Male 268 (52.4) 210 (52.0) 1.24 (0.47-3.30, p=0.665) 1.17 (0.44-3.15, p=0.752)
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.26 (0.90-1.76, p=0.182)
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.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)")
) %>%

# Tidy-up
select(-c(fit_id, index)) %>%
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) %>%
select(-fit_id, -index) %>%
dependent_label(apistrat, dependent) -> t
Dependent: API in 2000 (api00) Mean (sd) Coefficient (univariable) Coefficient (multivariable)
English language learners (percent)(ell) [0,84] 652.8 (121.0) -3.73 (-4.35–3.11, p<0.001) -0.48 (-1.25-0.29, p=0.222)
Meals eligible (percent)(meals) [0,100] 652.8 (121.0) -3.38 (-3.71–3.05, p<0.001) -3.14 (-3.70–2.59, p<0.001)
First year at the school (percent)(mobility) [1,99] 652.8 (121.0) -1.43 (-3.30-0.44, p=0.137) 0.23 (-0.54-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) %>%
select(-fit_id, -index) %>%
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.01, 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)