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
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
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
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
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 in chisq.test(age.factor, perfor.factor): Chi-squared approximation may
#> be incorrect
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
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
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
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 (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).
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 (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
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) |
|
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.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
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) |
|
Sex |
Female |
432 (47.9) |
13 (48.1) |
445 (47.9) |
1.000 |
|
Male |
470 (52.1) |
14 (51.9) |
484 (52.1) |
|
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
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 (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
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) |
|
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) |
|
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
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
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.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
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
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
|
|
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
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) |
|
Sex |
Female |
432 (47.9) |
13 (48.1) |
445 (47.9) |
1.000 |
|
Male |
470 (52.1) |
14 (51.9) |
484 (52.1) |
|
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
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) |
|
Sex |
Female |
432 (47.9) |
13 (48.1) |
445 (47.9) |
1.000 |
|
Male |
470 (52.1) |
14 (51.9) |
484 (52.1) |
|
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
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) |
|
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) |
|
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 levels_count
#> mort_5yr 2 "Alive", "Died", "(Missing)" 511, 404, 14
#> age.factor 3 "<40 years", "40-59 years", "60+ years" 70, 344, 515
#> sex.factor 2 "Female", "Male" 445, 484
#> obstruct.factor 2 "No", "Yes", "(Missing)" 732, 176, 21
#> levels_percent
#> mort_5yr 55.0, 43.5, 1.5
#> age.factor 7.5, 37.0, 55.4
#> sex.factor 48, 52
#> obstruct.factor 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.
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) |
|
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) |
|
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_explicit_na(mort_5yr)
) %>%
summary_factorlist(dependent, explanatory, p = TRUE, na_include = TRUE,
total_col = TRUE,
add_col_totals = TRUE, add_row_totals = TRUE) -> t
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) |
|
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) |
|
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
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) |
|
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) |
|
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.
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
#> Note: Using an external vector in selections is ambiguous.
#> ℹ Use `all_of(dependent)` instead of `dependent` to silence this message.
#> ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
#> Note: Using an external vector in selections is ambiguous.
#> ℹ Use `all_of(explanatory)` instead of `explanatory` to silence this message.
#> ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
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) |
|
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) |
|
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) |
|
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
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
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
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
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
Age (years) |
Mean (SD) |
29.7 (17.8) |
29.8 (17.5) |
Age |
<40 years |
33.7 (7.5) |
1.0 (7.4) |
|
40-59 years |
167.3 (37.3) |
4.7 (35.0) |
|
60+ years |
247.9 (55.2) |
7.8 (57.7) |
Sex |
Female |
216.7 (48.3) |
6.1 (44.7) |
|
Male |
232.2 (51.7) |
7.5 (55.3) |
Obstruction |
No |
350.9 (80.2) |
8.9 (65.3) |
|
Yes |
86.5 (19.8) |
4.7 (34.7) |
dim(colon_s) colon_s = colon_s %>% mutate(my_weights = runif(929, 0, 1))
Table 1 - Patient demographics —-
explanatory = c(“age”, “age.factor”, “sex.factor”, “obstruct.factor”) dependent = “perfor.factor”
colon_s %>% summary_factorlist(dependent, explanatory, p=TRUE)
colon_s %>% summary_factorlist(dependent, explanatory, p=TRUE, weights = “my_weights”, total_col = TRUE, digits = c(1, 1, 3, 1, 1))
colon_s %>% mutate(across(where(is.numeric) & any_of(explanatory), ~ .x * my_weights))
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
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
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
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
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
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
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
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
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
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
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
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
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
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...
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
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
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
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
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 %>%
lme4::glmer(mort_5yr ~ age.factor + (age.factor | hospital), family = "binomial", data = .) %>%
broom::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) %>%
na.omit() %>%
finalfit(dependent, explanatory) -> t
age.factor |
<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.factor |
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) |
obstruct.factor |
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) |
perfor.factor |
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
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
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-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–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-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-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-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
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-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–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-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-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-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
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
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
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
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
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
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
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) %>%
## Add univariable
ff_merge(
glmuni(colon_s, dependent, explanatory, weights = weights, family = quasibinomial) %>%
fit2df(estimate_suffix=" (univariable)")
) %>%
## Add multivariable
ff_merge(
glmmulti(colon_s, dependent, explanatory, weights = weights, family = quasibinomial) %>%
fit2df(estimate_suffix=" (multivariable)"),
last_merge = TRUE
) %>%
dependent_label(colon_s, dependent) -> t
Age |
<40 years |
31 (6.1) |
36 (8.9) |
- |
- |
|
40-59 years |
208 (40.7) |
131 (32.4) |
0.60 (0.34-1.05, p=0.073) |
0.65 (0.37-1.14, p=0.131) |
|
60+ years |
272 (53.2) |
237 (58.7) |
0.81 (0.47-1.38, p=0.439) |
0.89 (0.51-1.54, p=0.678) |
Sex |
Female |
243 (47.6) |
194 (48.0) |
- |
- |
|
Male |
268 (52.4) |
210 (52.0) |
1.03 (0.79-1.33, p=0.838) |
1.02 (0.78-1.34, p=0.870) |
Obstruction |
No |
408 (82.1) |
312 (78.6) |
- |
- |
|
Yes |
89 (17.9) |
85 (21.4) |
1.13 (0.81-1.57, p=0.468) |
1.15 (0.82-1.61, p=0.408) |
Perforation |
No |
497 (97.3) |
391 (96.8) |
- |
- |
|
Yes |
14 (2.7) |
13 (3.2) |
1.53 (0.70-3.39, p=0.284) |
1.46 (0.67-3.25, p=0.346) |
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
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
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
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
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
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
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) |