How to create a descriptive summary table (‘table 1’) using R

Intro

“Table 1”, that is a table providing the sample characteristics of an empirical study or clinical trial is an obligatory part of scientific publications. Since I started using {R} some ten years ago, I have come across a couple of packages and functions aiming to create such a table. In this blog post, I’ll give an overview about two of these packages and show how to use them.

Data

The data.frame I use is named ‘trial’ and is part of the {gtsummary} package. It contains the following variables:

head(gtsummary::trial)
## # A tibble: 6 x 8
##   trt      age marker stage grade response death ttdeath
##   <chr>  <dbl>  <dbl> <fct> <fct>    <int> <int>   <dbl>
## 1 Drug A    23  0.16  T1    II           0     0    24  
## 2 Drug B     9  1.11  T2    I            1     0    24  
## 3 Drug A    31  0.277 T1    II           0     0    24  
## 4 Drug A    NA  2.07  T3    III          1     1    17.6
## 5 Drug A    51  2.77  T4    III          1     1    16.4
## 6 Drug B    39  0.613 T4    I            0     1    15.6

Apparently, the data.frame contains a treatment variable (“trt”) with two categories (“Drug A” vs. “Drug B”) and several categorical and numerical variables.

{finalfit}

According to the package description, the {finalfit} package has the following purposes:

Generate regression results tables and plots in final format for publication. Explore models and export directly to PDF and ‘Word’ using ‘RMarkdown’.

However, with summary_factorlist(), the package also includes a function to create a table with summary statistics.

library(finalfit)
## Error in `library()`:
## ! there is no package called 'finalfit'
library(dplyr)
tab.ff <- gtsummary::trial %>%
  mutate(across(c(response, death),
    factor,
    levels = c(1, 0), labels = c("yes", "no")
  )) %>%
  summary_factorlist(
    dependent = "trt", # name of grouping / treatment variable
    explanatory = c("age", "", "marker", "stage", "response", "death", "ttdeath"),
    total_col = TRUE, # add column with statistics for the whole sample
    add_row_total = TRUE, # add column with number of valid cases
    include_row_missing_col = FALSE,
    na_include = TRUE # make variables' missing data explicit
  )
## Error in `summary_factorlist()`:
## ! could not find function "summary_factorlist"
tab.ff
## Error:
## ! object 'tab.ff' not found

Rather than printing a table, the summary_factorlist() function returns a data.frame which must be further processed and piped to a printing function. The following example shows how these steps may be done using the {labelled} and the {kableExtra} package.

library(labelled)
library(dplyr)
library(kableExtra, exclude = "group_rows")
gtsummary::trial %>%
  mutate(across(c(response, death),
    factor,
    levels = c(1, 0), labels = c("yes", "no")
  )) %>%
  # Add variable labels
  set_variable_labels(
    age = "Age [yrs]",
    marker = "Marker Level [ng/mL]",
    stage = "T Stage",
    grade = "Grade",
    response = "Tumor Response",
    death = "Patient Died",
    ttdeath = "Months to Death/Censor"
  ) %>%
  summary_factorlist(
    dependent = "trt", # name of grouping / treatment variable
    explanatory = c("age", "", "marker", "stage", "response", "death", "ttdeath"),
    total_col = TRUE, # add column with statistics for the whole sample
    add_row_total = TRUE, # add column with number of valid cases
    include_row_missing_col = FALSE,
    na_include = TRUE # make variables' missing data explicit
  ) %>%
  kbl(
    caption = "Baseline characteristics",
    booktabs = TRUE,
    col.names = c(
      " ", "Total N", " ",
      "Drug A", "Drug B", "Total"
    ),
    align = "lrlrrr",
  ) %>%
  kable_classic(full_width = FALSE)
## Error in `summary_factorlist()`:
## ! could not find function "summary_factorlist"

{gtsummary}

The description of the {gtsummary} package gives the following information:

The package creates presentation-ready tables summarizing data sets, regression models, and more. The code to create the tables is concise and highly customizable. Data frames can be summarized with any function, e.g. mean(), median(), even user-written functions. Regression models are summarized and include the reference rows for categorical variables. Common regression models, such as logistic regression and Cox proportional hazards regression, are automatically dentified and the tables are pre-filled with appropriate column headers.

For creating a table with summary statistics, the tbl_summary() function is required. In addition, the {dplyr} package should be loaded.

library(gtsummary)
library(dplyr)

The creation of the table works best using the pipe operator:

tbl.gts <- trial %>%
  # categorical variables must be factors or character strings
  mutate(across(c(response, death),
    factor,
    levels = c(1, 0), labels = c("yes", "no")
  )) %>%
  # apply the tbl_summary() function
  tbl_summary(
    by = trt, # Treatment variable
    label = list(
      age ~ "Age [yrs]",
      marker ~ "Marker Level [ng/mL]",
      stage ~ "T Stage",
      grade ~ "Grade",
      response ~ "Tumor Response",
      death ~ "Patient Died",
      ttdeath ~ "Months to Death/Censor"
    ),
    statistic = list(all_continuous() ~ "{mean} ({sd})"),
    # may be also ...
    # statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"),
    digits = all_continuous() ~ 1,
    missing_text = "(Missing)",
    include = everything() # select variables to be included into the table
  ) %>%
  add_overall() %>% # add column with statistics for the whole sample
  add_n() # add column with number of valid cases
tbl.gts

#tllrekyvog table {
font-family: system-ui, ‘Segoe UI’, Roboto, Helvetica, Arial, sans-serif, ‘Apple Color Emoji’, ‘Segoe UI Emoji’, ‘Segoe UI Symbol’, ‘Noto Color Emoji’;
-webkit-font-smoothing: antialiased;
-moz-osx-font-smoothing: grayscale;
}

#tllrekyvog thead, #tllrekyvog tbody, #tllrekyvog tfoot, #tllrekyvog tr, #tllrekyvog td, #tllrekyvog th {
border-style: none;
}

#tllrekyvog p {
margin: 0;
padding: 0;
}

#tllrekyvog .gt_table {
display: table;
border-collapse: collapse;
line-height: normal;
margin-left: auto;
margin-right: auto;
color: #333333;
font-size: 16px;
font-weight: normal;
font-style: normal;
background-color: #FFFFFF;
width: auto;
border-top-style: solid;
border-top-width: 2px;
border-top-color: #A8A8A8;
border-right-style: none;
border-right-width: 2px;
border-right-color: #D3D3D3;
border-bottom-style: solid;
border-bottom-width: 2px;
border-bottom-color: #A8A8A8;
border-left-style: none;
border-left-width: 2px;
border-left-color: #D3D3D3;
}

#tllrekyvog .gt_caption {
padding-top: 4px;
padding-bottom: 4px;
}

#tllrekyvog .gt_title {
color: #333333;
font-size: 125%;
font-weight: initial;
padding-top: 4px;
padding-bottom: 4px;
padding-left: 5px;
padding-right: 5px;
border-bottom-color: #FFFFFF;
border-bottom-width: 0;
}

#tllrekyvog .gt_subtitle {
color: #333333;
font-size: 85%;
font-weight: initial;
padding-top: 3px;
padding-bottom: 5px;
padding-left: 5px;
padding-right: 5px;
border-top-color: #FFFFFF;
border-top-width: 0;
}

#tllrekyvog .gt_heading {
background-color: #FFFFFF;
text-align: center;
border-bottom-color: #FFFFFF;
border-left-style: none;
border-left-width: 1px;
border-left-color: #D3D3D3;
border-right-style: none;
border-right-width: 1px;
border-right-color: #D3D3D3;
}

#tllrekyvog .gt_bottom_border {
border-bottom-style: solid;
border-bottom-width: 2px;
border-bottom-color: #D3D3D3;
}

#tllrekyvog .gt_col_headings {
border-top-style: solid;
border-top-width: 2px;
border-top-color: #D3D3D3;
border-bottom-style: solid;
border-bottom-width: 2px;
border-bottom-color: #D3D3D3;
border-left-style: none;
border-left-width: 1px;
border-left-color: #D3D3D3;
border-right-style: none;
border-right-width: 1px;
border-right-color: #D3D3D3;
}

#tllrekyvog .gt_col_heading {
color: #333333;
background-color: #FFFFFF;
font-size: 100%;
font-weight: normal;
text-transform: inherit;
border-left-style: none;
border-left-width: 1px;
border-left-color: #D3D3D3;
border-right-style: none;
border-right-width: 1px;
border-right-color: #D3D3D3;
vertical-align: bottom;
padding-top: 5px;
padding-bottom: 6px;
padding-left: 5px;
padding-right: 5px;
overflow-x: hidden;
}

#tllrekyvog .gt_column_spanner_outer {
color: #333333;
background-color: #FFFFFF;
font-size: 100%;
font-weight: normal;
text-transform: inherit;
padding-top: 0;
padding-bottom: 0;
padding-left: 4px;
padding-right: 4px;
}

#tllrekyvog .gt_column_spanner_outer:first-child {
padding-left: 0;
}

#tllrekyvog .gt_column_spanner_outer:last-child {
padding-right: 0;
}

#tllrekyvog .gt_column_spanner {
border-bottom-style: solid;
border-bottom-width: 2px;
border-bottom-color: #D3D3D3;
vertical-align: bottom;
padding-top: 5px;
padding-bottom: 5px;
overflow-x: hidden;
display: inline-block;
width: 100%;
}

#tllrekyvog .gt_spanner_row {
border-bottom-style: hidden;
}

#tllrekyvog .gt_group_heading {
padding-top: 8px;
padding-bottom: 8px;
padding-left: 5px;
padding-right: 5px;
color: #333333;
background-color: #FFFFFF;
font-size: 100%;
font-weight: initial;
text-transform: inherit;
border-top-style: solid;
border-top-width: 2px;
border-top-color: #D3D3D3;
border-bottom-style: solid;
border-bottom-width: 2px;
border-bottom-color: #D3D3D3;
border-left-style: none;
border-left-width: 1px;
border-left-color: #D3D3D3;
border-right-style: none;
border-right-width: 1px;
border-right-color: #D3D3D3;
vertical-align: middle;
text-align: left;
}

#tllrekyvog .gt_empty_group_heading {
padding: 0.5px;
color: #333333;
background-color: #FFFFFF;
font-size: 100%;
font-weight: initial;
border-top-style: solid;
border-top-width: 2px;
border-top-color: #D3D3D3;
border-bottom-style: solid;
border-bottom-width: 2px;
border-bottom-color: #D3D3D3;
vertical-align: middle;
}

#tllrekyvog .gt_from_md > :first-child {
margin-top: 0;
}

#tllrekyvog .gt_from_md > :last-child {
margin-bottom: 0;
}

#tllrekyvog .gt_row {
padding-top: 8px;
padding-bottom: 8px;
padding-left: 5px;
padding-right: 5px;
margin: 10px;
border-top-style: solid;
border-top-width: 1px;
border-top-color: #D3D3D3;
border-left-style: none;
border-left-width: 1px;
border-left-color: #D3D3D3;
border-right-style: none;
border-right-width: 1px;
border-right-color: #D3D3D3;
vertical-align: middle;
overflow-x: hidden;
}

#tllrekyvog .gt_stub {
color: #333333;
background-color: #FFFFFF;
font-size: 100%;
font-weight: initial;
text-transform: inherit;
border-right-style: solid;
border-right-width: 2px;
border-right-color: #D3D3D3;
padding-left: 5px;
padding-right: 5px;
}

#tllrekyvog .gt_stub_row_group {
color: #333333;
background-color: #FFFFFF;
font-size: 100%;
font-weight: initial;
text-transform: inherit;
border-right-style: solid;
border-right-width: 2px;
border-right-color: #D3D3D3;
padding-left: 5px;
padding-right: 5px;
vertical-align: top;
}

#tllrekyvog .gt_row_group_first td {
border-top-width: 2px;
}

#tllrekyvog .gt_row_group_first th {
border-top-width: 2px;
}

#tllrekyvog .gt_summary_row {
color: #333333;
background-color: #FFFFFF;
text-transform: inherit;
padding-top: 8px;
padding-bottom: 8px;
padding-left: 5px;
padding-right: 5px;
}

#tllrekyvog .gt_first_summary_row {
border-top-style: solid;
border-top-color: #D3D3D3;
}

#tllrekyvog .gt_first_summary_row.thick {
border-top-width: 2px;
}

#tllrekyvog .gt_last_summary_row {
padding-top: 8px;
padding-bottom: 8px;
padding-left: 5px;
padding-right: 5px;
border-bottom-style: solid;
border-bottom-width: 2px;
border-bottom-color: #D3D3D3;
}

#tllrekyvog .gt_grand_summary_row {
color: #333333;
background-color: #FFFFFF;
text-transform: inherit;
padding-top: 8px;
padding-bottom: 8px;
padding-left: 5px;
padding-right: 5px;
}

#tllrekyvog .gt_first_grand_summary_row {
padding-top: 8px;
padding-bottom: 8px;
padding-left: 5px;
padding-right: 5px;
border-top-style: double;
border-top-width: 6px;
border-top-color: #D3D3D3;
}

#tllrekyvog .gt_last_grand_summary_row_top {
padding-top: 8px;
padding-bottom: 8px;
padding-left: 5px;
padding-right: 5px;
border-bottom-style: double;
border-bottom-width: 6px;
border-bottom-color: #D3D3D3;
}

#tllrekyvog .gt_striped {
background-color: rgba(128, 128, 128, 0.05);
}

#tllrekyvog .gt_table_body {
border-top-style: solid;
border-top-width: 2px;
border-top-color: #D3D3D3;
border-bottom-style: solid;
border-bottom-width: 2px;
border-bottom-color: #D3D3D3;
}

#tllrekyvog .gt_footnotes {
color: #333333;
background-color: #FFFFFF;
border-bottom-style: none;
border-bottom-width: 2px;
border-bottom-color: #D3D3D3;
border-left-style: none;
border-left-width: 2px;
border-left-color: #D3D3D3;
border-right-style: none;
border-right-width: 2px;
border-right-color: #D3D3D3;
}

#tllrekyvog .gt_footnote {
margin: 0px;
font-size: 90%;
padding-top: 4px;
padding-bottom: 4px;
padding-left: 5px;
padding-right: 5px;
}

#tllrekyvog .gt_sourcenotes {
color: #333333;
background-color: #FFFFFF;
border-bottom-style: none;
border-bottom-width: 2px;
border-bottom-color: #D3D3D3;
border-left-style: none;
border-left-width: 2px;
border-left-color: #D3D3D3;
border-right-style: none;
border-right-width: 2px;
border-right-color: #D3D3D3;
}

#tllrekyvog .gt_sourcenote {
font-size: 90%;
padding-top: 4px;
padding-bottom: 4px;
padding-left: 5px;
padding-right: 5px;
}

#tllrekyvog .gt_left {
text-align: left;
}

#tllrekyvog .gt_center {
text-align: center;
}

#tllrekyvog .gt_right {
text-align: right;
font-variant-numeric: tabular-nums;
}

#tllrekyvog .gt_font_normal {
font-weight: normal;
}

#tllrekyvog .gt_font_bold {
font-weight: bold;
}

#tllrekyvog .gt_font_italic {
font-style: italic;
}

#tllrekyvog .gt_super {
font-size: 65%;
}

#tllrekyvog .gt_footnote_marks {
font-size: 75%;
vertical-align: 0.4em;
position: initial;
}

#tllrekyvog .gt_asterisk {
font-size: 100%;
vertical-align: 0;
}

#tllrekyvog .gt_indent_1 {
text-indent: 5px;
}

#tllrekyvog .gt_indent_2 {
text-indent: 10px;
}

#tllrekyvog .gt_indent_3 {
text-indent: 15px;
}

#tllrekyvog .gt_indent_4 {
text-indent: 20px;
}

#tllrekyvog .gt_indent_5 {
text-indent: 25px;
}

#tllrekyvog .katex-display {
display: inline-flex !important;
margin-bottom: 0.75em !important;
}

#tllrekyvog div.Reactable > div.rt-table > div.rt-thead > div.rt-tr.rt-tr-group-header > div.rt-th-group:after {
height: 0px !important;
}

Characteristic N Overall
N = 200
1
Drug A
N = 98
1
Drug B
N = 102
1
Age [yrs] 189 47.2 (14.3) 47.0 (14.7) 47.4 (14.0)
    (Missing) 11 7 4
Marker Level [ng/mL] 190 0.9 (0.9) 1.0 (0.9) 0.8 (0.8)
    (Missing) 10 6 4
T Stage 200
    T1 53 (27%) 28 (29%) 25 (25%)
    T2 54 (27%) 25 (26%) 29 (28%)
    T3 43 (22%) 22 (22%) 21 (21%)
    T4 50 (25%) 23 (23%) 27 (26%)
Grade 200
    I 68 (34%) 35 (36%) 33 (32%)
    II 68 (34%) 32 (33%) 36 (35%)
    III 64 (32%) 31 (32%) 33 (32%)
Tumor Response 193 61 (32%) 28 (29%) 33 (34%)
    (Missing) 7 3 4
Patient Died 200 112 (56%) 52 (53%) 60 (59%)
Months to Death/Censor 200 19.6 (5.3) 20.2 (5.0) 19.0 (5.5)
1 Mean (SD); n (%)


Since the {gtsummary} packages contains functions to convert the {gtsummary} object to object types required by other popular table-specific R packages, e.g. as_kable_extra(), as_flextable() etc., {gtsummary} tables can be easily rendered as .docx, .pdf or .html. The following example shows how to print a {gtsummary} object using the {kableExtra} package.

library(kableExtra)
tbl.gts %>%
  as_kable_extra(
    caption = "Baseline characteristics",
    booktabs = TRUE,
    align = "lcccc",
  ) %>%
  kable_classic(full_width = FALSE)
Baseline characteristics
Characteristic N Overall
N = 200
Drug A
N = 98
Drug B
N = 102
Age [yrs] 189 47.2 (14.3) 47.0 (14.7) 47.4 (14.0)
(Missing) 11 7 4
Marker Level [ng/mL] 190 0.9 (0.9) 1.0 (0.9) 0.8 (0.8)
(Missing) 10 6 4
T Stage 200
T1 53 (27%) 28 (29%) 25 (25%)
T2 54 (27%) 25 (26%) 29 (28%)
T3 43 (22%) 22 (22%) 21 (21%)
T4 50 (25%) 23 (23%) 27 (26%)
Grade 200
I 68 (34%) 35 (36%) 33 (32%)
II 68 (34%) 32 (33%) 36 (35%)
III 64 (32%) 31 (32%) 33 (32%)
Tumor Response 193 61 (32%) 28 (29%) 33 (34%)
(Missing) 7 3 4
Patient Died 200 112 (56%) 52 (53%) 60 (59%)
Months to Death/Censor 200 19.6 (5.3) 20.2 (5.0) 19.0 (5.5)
1 Mean (SD); n (%)

Other packages

Other packages I have used in the past are {qwraps2} and {tableone}. However, the {gtsummary} package seems to offer the most convenient way to create a “table one”.

Statistische Grundlagen mit R: Streuungsmaße

Nachdem ich mich im ersten Teil meiner Serie “Statistische Grundlagen mit R” mit Maßen der zentralen Tendenz beschäftigt habe, erkläre ich heute die Bedeutung verschiedener Steuungsmaße.

Streuungsmaße

Um einen Datenvektor möglichst treffend zu beschreiben, reicht die Kenntnis der Zentralitätsmaße in der Regel nicht aus. Stellen wir uns z.B. vor, dass an einer Grundschule ein Altpapiersammelwettbewerb durchgeführt wird, bei dem die drei Schulklassen, die am meisten Papier gesammelt haben Preise bekommen: die erstplatzierte Klasse 100€, die zweitplatzierte Klasse 50€ und die drittplatzierte Klasse 20€. Wenn wir uns vorstellen, dass die erstplatzierte Klasse insgesamt 800 kg, die zweitplatzierte Klasse 700 kg und die drittplatzierte Klasse 600 kg gesammelt habt, würde uns die Vergabe der Preise im Großen und Ganzen gerecht erscheinen. Stellen wir uns hingegen vor, dass die erstplatzierte Klasse 701 kg, die zweitplatzierte Klasse genau 700 kg und die drittplatzierte Klasse 699 kg Papier gesammelt hat, erscheint die Preisvergabe weniger angemessen zu sein, weil die Unterschiede zwischen den gesammelte Altpapiermengen äußerst klein sind und möglicherweise sogar auf Messungenauigkeiten (beim Wiegen des Papiers) zurückzuführen sind. Berechnen wir den Mittelwert für beide Szenarien, sehen wir jedoch keinen Unterschied; beide Mittelwerte sind identisch:

mean(c(600, 700, 800)) == mean(c(699, 700, 701))
## [1] TRUE

Mit Streuungsmaßen lässt sich, wie der Name schon sagt, die Streuung bzw. Entfernung der einzelnen Datenpunkte von einem Zentralitätsmaß beschreiben.

Spanne

Das wohl einfachste Streuungsmaß ist die Spanne (engl.: Range). Diese lässt sich berechnen, indem man den kleinsten Wert vom größen des Datenvektors abzieht:

spanne1 <- 800 - 600
spanne2 <- 701 - 699

Obwohl beide Datenvektoren denselben Mittelwert haben, unterscheiden sich die Spannen deutlich (200 vs. 2). Die Spanne ist als Streuungsmaß allerdngs nur sehr eingeschränkt zu gebrauchen, weil sie sehr stark durch Extremwerte beeinflusst wird. So könnte man sich vorstellen, dass bei dem besagten Altpapiersammelwettbewerb 8 von 10 Schulklassen zwischen 400 und 500 kg Altpapier sammeln, eine Klasse 1.000 kg und eine Klasse überhaupt kein Papier. Obwohl die meisten Datenpunkte zwischen 400 und 500 kg liegen, nimmt die Spanne den Wert 1.000 (1.000 – 0) an.

Abweichungen vom Mittelwert

Eine weitere Möglichkeit zur Beschreibung der Streuung eines Datenvektors ist die Abweichung der einzelnen Datenpunkte vom Mittelwert aller Datenpunkte. Bezugnehmend auf unser Beispiel, weichen im ersten Szenario zwei von drei Datenpunkten um jeweils 100 kg vom Mittelwert ab, während der dritte Wert mit dem Mittelpunkt identisch ist. Bildet man die Summe dieser Abweichungen, erhalten wir das folgende Ergebnis:

-100 + 0 + 100
## [1] 0

Berechnet man die Streuung nach dieser Gleichung, muss diese zwangsläufig den Wert 0 ergeben, da sich positive und negative Differenzen gegenseitig aufheben. Um dies zu verhindern, ließen sich z.B. die absoluten Beträge der Differenzen aufaddieren. In der Statistik entledigt man sich des Minuszeichens jedoch durch eine Quadrierung der einzelnen Differenzen.

Summe der quatrierten Abweichungen

Die einzelnen quatrierten Differenzen werden zu einer Quadratsumme zusammengefasst. Für unsere beiden Szenarien der Altpapiersammelaktion lassen sich die Quadratsummen folgendermaßen berechnen:

(qs1 <- (-100)^2 + 0^2 + 100^2)
## [1] 20000
(qs2 <- (-1)^2 + 0^2 + 1^2)
## [1] 2

Als Streuungsmaß hat die Quadratsumme jedoch den Nachteil, dass ihr Betrag mit steigender Anzahl an Datenpunkten zunimmt.

Varianz

Um eine Abhängikeit des Streungsmaßes von der Anzahl der Datenpunkte zu verhindern, lässt sich die Summe der quatrierten Abweichungen durch die Anzahl der Datenpunkte (Beobachtungen) teilen. Mit dieser Gleichung haben wir die Varianz eines Datenvektors definiert. Um zu demonstrieren, wie sich die Varianz eines Datenvektors mit R berechnen lässt, kommen wir wieder auf unser Beispiel vom Altpapiersammelwettbewerb zurück:

altpapier <- c(600, 700, 800) # vgl. Szenario 1
(altpapier - mean(altpapier)) ^ 2 # Liste quatrierter Abweichungen vom Mittelwert
## [1] 10000     0 10000
sum((altpapier - mean(altpapier)) ^ 2) # Summe der quatrierten Abweichungen
## [1] 20000
(varianz <- sum((altpapier - mean(altpapier)) ^ 2) / length(altpapier)) # Teile durch Anzahl der Beobachtungen
## [1] 6666.667

Standardabweichung

Die Interpretation der Varianz gestaltet sich jedoch schwierig, da sie die Einheiten des Datenvektors (in unserem Beispiel Kilogramm) quatriert widergibt. So hat der Vektor altpapier eine Varianz von etwa 6666.7 kg². Um die Standardabweichung zu berechnen, braucht man nun nur noch die Quatratwurzel der Varianz ziehen:

sqrt(varianz)
## [1] 81.64966

Im Gegensatz zur Varianz, hat die Standardabweichung dieselbe Einheit wie der Datenvektor, in unserem Fall etwa 81.6 kg (Szenario 1) und 0.816 kg (Szenario 2).

Empfohlene Literatur

  • Jeffrey M. Stanton, Reasoning with Data: An Introduction to Traditional and Bayesian Statistics Using R, Guilford Press: 2017.

Statistische Grundlagen mit R: Maße der zentralen Tendenz

Beginnend mit dem heutigen Tag veröffentliche ich auf “Scripts and Statistics” einige Beiträge unter der Überschrift “Statistische Grundlagen mit R”. Im ersten Teil erkläre ich die drei wichtigsten Maße der zentralen Tendenz.

Maße der zentralen Tendenz

Mit Maßen der zentralen Tendenz lässt sich ein Datenvektor (eine Reihe von Zahlen) mit einem Kennwert charakterisieren. Die drei häufigsten Zentralitätsmaße sind Mittelwert, Median und Modus.

Mittelwert

Das wohl bekannteste statistische Maß ist das arithmetische Mittel, auch Mittelwert genannt (engl.: “mean”). Der Mittelwert wird berechnet, indem man alle Werte eines Datenvektors addiert und die Summe durch die Anzahl der Einzelbeobachtungen (Länge des Vektors) teilt. So lässt sich z.B. die Anzahl der an einem Spieltag der Fußball-Bundesliga erzielten Treffer pro Spiel wie folgt darstellen:

treffer <- c(2, 3, 6, 0, 1, 3, 9, 3, 2)

Mit der c()-Funktion wird ein Datenvektor erzeugt und unter der Bezeichnung treffer abgespeichert. Um den Mittelwert zu berechnen, müssen alle Zahlen dieses Vektors addiert und durch 9 (Anzahl der Spiele pro Spieltag) dividiert werden:

sum(treffer) / length(treffer) # Ausführliche Berechnung
## [1] 3.222222
mean(treffer) # Mit R-Funktion berechnet
## [1] 3.222222

In der ersten Zeile wird der Mittelwert – wie soeben erklärt – berechnet, in der zweiten Zeile mit der mean()-Funktion. Beide Rechenwege kommen auf dasselbe Ergebnis.

Median

Der Median (engl.: “median”) lässt sich etwas salopp auch als “Halbe-Strecke-Wert” (halfway value) bezeichnen. Ordnet man einen Datenvektor aufsteigend vom kleinsten bis zum größten Wert und sucht sich den Wert, der sich genau in der Mitte befindet, so erhält man den Median. Mit der sort()-Funktion wird der Datenvektor aufsteigend geordnet:

sort(treffer)
## [1] 0 1 2 2 3 3 3 6 9

Da der Datenvektor treffer aus neun Einzelwerten besteht, muss der fünfte und damit mittlere Wert dem Median entsprechen:

sort(treffer)[5] # Median als mittlerer Wert des sortierten Datenvektors
## [1] 3
median(treffer) # Mit R-Funktion berechnet
## [1] 3

Beide Werte sind identisch. Wie wir sehen, ist der Median mit 3 Treffern ewtas kleiner als der Mittelwert mit 3.2222222 Treffern. Dieses Ergebnis weist auf einen wichtigen Vorteil des Medians hin: dieser ist robuster als der Mittelwert in Hinblick auf Extremwerte (engl.: “outliers”). Ein Extremwert ist ein Wert, der deutlich höher oder niedriger ist als die meisten anderen Werte. In unserem Beispiel sind die 9 Tore, die bei einem Bundesligaspiel erzielt wurden, ein Extremwert.

Modus

Der Modus (engl.: “mode”) ist der häufigste Wert des Datenvektors. Tritt ein bestimmter Wert eines Datenvektors besonders häufig auf, so lässt sich mit dem Modus der typischste Wert dieses Vektors angeben. Der Modus ist in noch stärkerem Maße als der Median robust gegenüber Extremwerten. Da in den R-Basispaketen keine Funktion zur Berechnung des Modus implementiert ist, muss ein zusätzliches R-Paket (DescTools) installiert und geladen werden. Mit der Mode()-Funktion dieses Pakets lässt sich nun der Modus unseres Datenvektors berechnen:

install.packages('DescTools')
library(DescTools)
Mode(treffer)
## [1] 3

Der Modus zeigt an, dass pro Bundesligaspiel am häufigsten 3 Tore fallen.

Resumee

Anhand der soeben erläuterten Maße der zentralen Tendenz ist eine zusammenfassende Beschreibung eines Datenvektors möglich. Der Mittelwert liefert uns einen präzisen, arithmetischen Mittelwert des Datenvektors, ist aber vergleichsweise anfällig gegenüber Extremwerten. Den Median kann man sich gut als mittleren Punkt einer Strecke vorstellen, der diese Strecke in zwei gleich große Hälften teilt. Wenn sich Mittelwert und Median unterscheiden, ist dies ein Hinweis darauf, dass der Datenvektor Extremwerte beinhaltet, die sich von den gewöhnlichen Werten stark unterscheiden. Der Modus zeigt uns, welcher Wert des Datenvektors am häufigsten vorkommt.

Empfohlene Literatur

  • Jeffrey M. Stanton, Reasoning with Data: An Introduction to Traditional and Bayesian Statistics Using R, Guilford Press: 2017.

Scoring the PHQ-9 Questionnaire Using R

Intro

The PHQ-9 is the nine-item depression module of the Patient Health Questionnaire. Each of the items is scored on a 4-point Likert scale ranging from 0 (not at all) to 3 (nearly every day). The items are summed to obtain a total score ranging from 0 to 27 with higher scores indicating greater severity of depression. Based on the total score, different levels of severity can be evaluated with 0–4, 5–9, 10–14, 15–19 and 20–27 points indicating “minimal”, “mild”, “moderate”, “moderately severe” and “severe” depression.

The PHQ-9 questionnaire may be found under the following link.

In this blog post, I show how to calculate the PHQ-9 score and the PHQ-9 severety levels.

Packages and data

The dataset we are going to use was published in Plos One. The file has got a Digital Object Identifier (doi) and may be imported into R using the read_delim() function of the readr package.

library(readr)
library(dplyr)
library(ggplot2)

df.phq9 <- readr::read_delim("https://doi.org/10.1371/journal.pone.0156167.s001", 
                             delim = ";", 
                             escape_double = FALSE, 
                             trim_ws = TRUE) %>%
            select(starts_with('phq9'))

glimpse(df.phq9)
## Observations: 1,337
## Variables: 9
## $ phq9_1 <int> 1, 3, 2, 0, 0, 0, 1, 0, 0, 2, 1, 1, 0, 3, 0, 0, 0, 2, 0...
## $ phq9_2 <int> 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0...
## $ phq9_3 <int> 3, 2, 2, 2, 1, 0, 1, 3, 1, 0, 1, 1, 0, 3, 1, 0, 0, 0, 0...
## $ phq9_4 <int> 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 1, 3, 0, 1, 0, 0, 0, 1, 0...
## $ phq9_5 <int> 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 2, 0, 1, 0, 0, 0, 0, 0...
## $ phq9_6 <int> 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0...
## $ phq9_7 <int> 0, 1, 1, 1, 0, 1, 0, 0, 0, 3, 1, 1, 0, 1, 0, 0, 0, 0, 0...
## $ phq9_8 <int> 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0...
## $ phq9_9 <int> 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

The Scoring Function

The scoring_phq9 function requires a data frame containing the PHQ-9 items (data) and a vector containing the items' names (items.phq9) as input parameters.

scoring_phq9 <- function(data, items.phq9) {
  data %>%
    mutate(nvalid.phq9 = rowSums(!is.na(select(., items.phq9))),
           nvalid.phq9 = as.integer(nvalid.phq9),
           mean.temp = rowSums(select(., items.phq9), na.rm = TRUE)/nvalid.phq9,
           phq.01.temp = as.integer(unlist(data[items.phq9[1]])),
           phq.02.temp = as.integer(unlist(data[items.phq9[2]])),
           phq.03.temp = as.integer(unlist(data[items.phq9[3]])),
           phq.04.temp = as.integer(unlist(data[items.phq9[4]])),
           phq.05.temp = as.integer(unlist(data[items.phq9[5]])),
           phq.06.temp = as.integer(unlist(data[items.phq9[6]])),
           phq.07.temp = as.integer(unlist(data[items.phq9[7]])),
           phq.08.temp = as.integer(unlist(data[items.phq9[8]])),
           phq.09.temp = as.integer(unlist(data[items.phq9[9]]))) %>%
    mutate_at(vars(phq.01.temp:phq.09.temp),
              funs(ifelse(is.na(.), round(mean.temp), .))) %>%
    mutate(score.temp = rowSums(select(., phq.01.temp:phq.09.temp), na.rm = TRUE),
           score.phq9 = ifelse(nvalid.phq9 >= 7, as.integer(round(score.temp)), NA),
           cutoff.phq9 = case_when(
             score.phq9 >= 20 ~ 'severe',
             score.phq9 >= 15 ~ 'moderately severe',
             score.phq9 >= 10 ~ 'moderate',
             score.phq9 >= 5 ~ 'mild',
             score.phq9 < 5 ~ 'minimal'),
             cutoff.phq9 = factor(cutoff.phq9, levels = c('minimal', 'mild',
                                                          'moderate', 'moderately severe',
                                                          'severe'))) %>%
    select(-ends_with("temp"))

}

Example

The function adds three variables to the original data frame:

  • nvalid.phq9: Number of variables with valid values,
  • score.phq9: PHQ-9 score (0 – 27),
  • cutoff.phq9: PHQ-9 severety levels (minimal, mild, moderate, moderately severe, severe)
items.phq9 <- paste0('phq9_', seq(1, 9, 1))
df.phq9 <- df.phq9 %>%
  scoring_phq9(., items.phq9)
glimpse(df.phq9)
## Observations: 1,337
## Variables: 12
## $ phq9_1      <int> 1, 3, 2, 0, 0, 0, 1, 0, 0, 2, 1, 1, 0, 3, 0, 0, 0,...
## $ phq9_2      <int> 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,...
## $ phq9_3      <int> 3, 2, 2, 2, 1, 0, 1, 3, 1, 0, 1, 1, 0, 3, 1, 0, 0,...
## $ phq9_4      <int> 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 1, 3, 0, 1, 0, 0, 0,...
## $ phq9_5      <int> 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 2, 0, 1, 0, 0, 0,...
## $ phq9_6      <int> 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,...
## $ phq9_7      <int> 0, 1, 1, 1, 0, 1, 0, 0, 0, 3, 1, 1, 0, 1, 0, 0, 0,...
## $ phq9_8      <int> 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0,...
## $ phq9_9      <int> 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ nvalid.phq9 <int> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,...
## $ score.phq9  <int> 7, 10, 7, 9, 3, 2, 3, 4, 5, 7, 7, 8, 0, 11, 1, 0, ...
## $ cutoff.phq9 <fct> mild, moderate, mild, mild, minimal, minimal, mini...

Visualization

PHQ-9 Score

ggplot(df.phq9, aes(score.phq9)) +
  geom_density(fill = 'blue', alpha = 0.2) +
  scale_x_continuous(limits = c(0, 27), breaks = c(0,5,10,15,20,27)) +
  labs(x = 'PHQ-9 Score', y = 'Density') +
  theme_bw()

plot of chunk unnamed-chunk-4

PHQ-9 Severety Levels

ggplot(df.phq9, aes(x = cutoff.phq9, fill = cutoff.phq9)) +
  geom_bar(colour = 'black') +
  scale_fill_brewer(type = 'seq') +
  labs(x = NULL, y = NULL, fill = NULL) +
  theme_bw()

plot of chunk unnamed-chunk-5

Sample Size Calculation Using R

Preample

“Sample Size Calculation Using R” – the title of this blog post – sounds very comprehensive. However, this post only deals with the sample size collection for one particular test (two sample t-test). A collection of functions for basic sample size calculation can be found in the pwr package.

Introduction

A couple of months ago, I attended some statistics course about sample size calculation and sample size adjustment. Amongst others, we spoke about a randomized placebo-controlled trial with patients suffering of depression. The trial's primary end point was the difference between placebo and treatment group in the HAM-D total score (some depression score) between baseline and the end of therapy. During the course, the professor asked us to write a “two-liner” returning the sample size per group for a two-sided, two sample t-test with the following parameters:

  • alpha = 0.05,
  • sigma = 8 ,
  • delta = 4 and
  • power ≤ 0.90

Here is my solution:

In a first step, I wrote a loop returning the power rising from 0.1 to 0.9 and the corresponding sample sizes respectively. The loop may either be a repeat or a while loop.

While Loop

A while loop runs and repeats while a specified condition returns TRUE. The loop terminates only when the condition evaluates to FALSE. In our case, the loop stops evaluating as soon as p > 0.9.

p <- 0.05
while(p <= 0.9){
  p <- p + 0.05
  t <- power.t.test(sd = 8, delta = 4, sig.level=0.05, power = p, type="two.sample", alternative = "two.sided")
  print(paste(t$power, ceiling(t$n), sep = ' - '))
}

Repeat Loop

In a repeat loop it is required to define a so-called break declaration to stop repeating the code. Just like the while loop, our repeat loop stops evaluating as soon as p > 0.9.

Both, the while and the repeat loop return the same result.

p <- 0.05
repeat{
  p <- p + 0.05
  t <- power.t.test(sd = 8, delta = 4, sig.level=0.05, power = p, type="two.sample", alternative = "two.sided")
  print(paste(t$power, ceiling(t$n), sep = ' - '))
  if(p > 0.90){
    break
  }
}
## [1] "0.1 - 5"
## [1] "0.15 - 8"
## [1] "0.2 - 12"
## [1] "0.25 - 15"
## [1] "0.3 - 18"
## [1] "0.35 - 21"
## [1] "0.4 - 25"
## [1] "0.45 - 28"
## [1] "0.5 - 32"
## [1] "0.55 - 36"
## [1] "0.6 - 41"
## [1] "0.65 - 45"
## [1] "0.7 - 51"
## [1] "0.75 - 57"
## [1] "0.8 - 64"
## [1] "0.85 - 73"
## [1] "0.9 - 86"

For a very good introduction to loops, see Davies (2016): The Book of R. No Starch Press: San Francisco.

Storing the Results

In case we want to use the results for whatever purpose, we need to store them in a vector. The following example heavily borrows from StackOverflow. The loop returns a data frame containing two variables: pwr (power) and num (sample size).

p <- 0.05
pwr <- c()
num <- c()
while(p <= 0.9){
  p <- p + 0.05
  t <- power.t.test(sd = 8, delta = 4, sig.level=0.05, power = p, type="two.sample", alternative = "two.sided")
  pwr <- c(pwr, t$power)
  num <- c(num, ceiling(t$n))
  df <- data.frame(pwr = pwr, n = num)
}
rm(list = c(setdiff(ls(), c("df"))))
dplyr::glimpse(df)
## Observations: 17
## Variables: 2
## $ pwr <dbl> 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55...
## $ n   <dbl> 5, 8, 12, 15, 18, 21, 25, 28, 32, 36, 41, 45, 51, 57, 64, ...

Plotting Power and Sample Size

Finally, we plot the relation between power and sample size using the ggplot2 package. The theme I use, stems from the ggpubr package.

ggplot(df, aes(pwr, n)) +
  geom_point() +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.1)) +
  scale_y_continuous(breaks = seq(0, max(df$n)+10, 10)) +
  geom_text(aes(label = n), nudge_y = 4) +
  ggpubr::theme_pubr() +
  labs(x = 'Power', y = 'Sample Size', 
       title = 'Power and Sample Size', 
        subtitle = expression('Two-sided, two sample t-test with' ~ sigma ~ '= 8,' ~ delta ~ '= 4,' ~ alpha ~ '= 0.05'))

plot of chunk unnamed-chunk-4

As we see, a power of 0.90 requires a sample size of 86 patients per treatment group for a two-sided, two sample t-test with sd = 8, delta = 4, and alpha = 0.05.

Multiple Imputations of missing values — ein Erfahrungsbericht

Intro

In der Forschungspraxis werden Datenpunkte mit fehlenden Werten oft von der Analyse ausgeschlossen oder durch einen zentralen Lageparameter (Mittelwert, Median) ersetzt. Als Alternative dazu hat sich das Verfahren der Multiplen Imputation etabliert, welches in diesem Beitrag erklärt wird.

Die im Folgenden geschilderten Erfahrungen habe ich mit der Statistiksoftware R und dem mice-Paket gesammelt [2].

R Pakete

In dieser Blog-Post habe ich die folgenden R-Pakete verwendet:

library(dplyr)
library(ggplot2)
library(mice)

Hintergrund

Bei der Durchführung statistischer Auswertungen ist man oft mit dem Auftreten fehlender Werte konfrontiert. Schließt man Datenpunkte mit fehlenden Werten von der Analyse aus, führt das zu einer Reduktion der Fallzahl und damit zu einer geringeren Power der statistischen Tests. Da Werte nur selten zufällig fehlen, kann ein Fallausschluss zu einer Verzerrung der Ergebnisse führen. Ersetzt man fehlende Werte hingegen durch den Mittelwert oder Median der entsprechenden Messgröße führt dies zu einer Reduktion der Varianz, womit sich das Risiko für einen alpha-Fehler erhöht.

Als Alternative zu diesen und anderen Methoden wurde in den letzten Jahrzehnten die Methode der Multiplen Imputation (MI) entwickelt [5]. Ausgehend von einem Datensatz mit fehlenden Werten wird eine festgelegte Anzahl von kompletten Datensätzen (m) erzeugt, indem die fehlenden Werte durch sogenannte plausible Werte ersetzt werden. Diese plausiblen Werte werden einer Datenverteilung entnommen, die für die jeden fehlenden Wert modelliert wird. Die m imputierten Datensätze gleichen sich hinsichtlich der beobachteten Datenpunkte, unterscheiden sich jedoch in den imputierten fehlenden Werten, wobei die Unterschiede der Unsicherheit des Imputationsprozesses Rechnung tragen [1]. Das Verfahren der MI lässt sich mittlerweile mit (fast) allen gängigen Softwarepaketen (R, SAS, SPSS, STATA) durchführen.

Vorgehen

Entscheidend für die Qualität einer MI ist die Spezifikation des Imputationsmodells, welche in Abhängigkeit vom Variablentyp erfolgen muss. Für die Imputation numerischer Variablen kann z.B. die Methode des “Predictive mean matching” (pmm) verwendet werden. Diese Methode hat den Vorteil, dass die imputierten Werte nur innerhalb der durch die Variable vorgegebenen Spannweite liegen. Imputiert man z.B. die fehlenden Werte der Skala eines EORTC-Fragebogens, liegen die imputierten Werte im Bereich 0 bis 100.

Datensätze, in denen Längsschnittdaten erfasst werden, können in zwei verschiedenen Datenformaten vorliegen. So ist in einem flachen Datensatz für jeden Fall (z.B. Patienten) genau eine Zeile reserviert, während pro Messgröße (z.B. Lebensqualitätsskala) und Befragungszeitpunkt jeweils eine Variable existiert. In einem tiefen Datensatz hingegen entspricht die Anzahl der Zeilen pro Fall der Anzahl der Messzeitpunkte, während pro Messgröße nur eine Variable angelegt wird, welche die Messungen aller Messzeitpunkte enthält. Der Zeitpunkt der Intervention bzw. Befragung wird in einer separaten Variable angegeben. Für die Durchführung einer MI empfiehlt es sich, mit einem tiefen Datensatz zu arbeiten.

Arbeitet man mit Lebensqualitätsfragebögen, so steht man vor der Entscheidung, die Imputationen entweder auf Ebene der Einzelitems oder daraus berechneter Skalen durchzuführen. Da statistische Analysen in der Regel auf der Grundlage von Skalen erfolgen, empfiehlt es sich, die fehlenden Skalenwerte zu imputieren. Dazu kommt, dass die Anzahl der Skalen deutlich geringer ist als die Anzahl der Einzelitems. Mit der Anzahl der zu imputierenden Variablen reduziert sich auch die für die MI benötigte Zeit bzw. Rechnerleistung [4]. Darüber hinaus erscheint es sinnvoll, alle Variablen mit fehlenden Werten in einem “Durchgang” zu imputieren.

Probleme

Imputiert man fehlende Werte eines Datensatzes in tiefer Form, erhält man ein Objekt der Klasse mids. Dieses enthält neben den ursprünglichen, nicht-imputierten Daten auch die Daten der m imputierten Datensätze. Grundlage für die Datenanalyse ist eine gepoolte Version der imputierten Datensätze. Benötigt man aber für die Datenanalyse den Baselinewert einer Variable als eigene Variable, so steht man vor dem Problem, die Daten umformen zu müssen. Da dies unter Beibehaltung der m imputierten Datensätze geschehen soll, muss dass mids-Objekt zunächst in eine gestaplete Matrix (auch stacked oder tall matrix) umgewandelt werden [3]. Hat man diese Matrix in die gewünschte Form gebracht, so kann sie problemlos in ein mids-Objekt zurückverwandelt werden, auf dessen Grundlage statistische Modelle berechnet werden können.

Eigene Erfahrungen haben außerdem gezeigt, dass eine zu starke Collinearität der imputierten Variablen zu einer Verzerrung der imputierten Werte führen kann. Möchte man beispielsweise fehlende Werte von Sub-Skalen eines Fragebogens sowie der entsprechenden Gesamtskala imputieren, empfielt es sich, die Gesamtskala als informations-gebende Variable aus dem Imputationsprozess auszuschließen.

Güteüberprüfung

Um die Güte der MI beurteilen zu können sollte geprüft werden, ob die Verteilung der nicht-imputierten in etwa mit jener der imputierten Werte übereinstimmt. Dies kann z.B. mit Hilfe von Histogrammen oder einer Kombination aus Stripcharts und Violinplots erfolgen* (vgl. Abb. 1).

Abb. 1: Stripchart zur Güteüberprüfung der MI

Abb. 1 zeigt eine Kombination aus Stripchart und Violinplot für eine Variable (n=200) mit 50 fehlenden Werten, die 5 mal imputiert wurden (Methode: pmm). Es wird deutlich, dass die Verteilung der gepoolten imputierten Werte in etwa der Verteilung der nicht-imputierten Werte entspricht.

Fazit für die Praxis

Im Vergleich zum traditionellen Umgang mit fehlenden Werten (Fallausschluss, einfache Imputationen) führen Multiple Imputationen weder zu einem Verlust statistischer Power noch zu einer Unterschätzung der Varianz. Auch die Schätzwerte werden in der Regel nicht zu stark verzerrt. In der Forschungspraxis bieten Multiplen Imputationen eine gute Möglichkeit im Umgang mit fehlenden Werten.

Literatur

  1. Buuren, Stef van. 2012. Flexible Imputation of Missing Data. CRC Press.

  2. Buuren, Stef, and Karin Groothuis-Oudshoorn. 2011. Mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software 45 (3). http://doc.utwente.nl/78938/.

  3. Errickson, Josh. 2017. Multiple Imputation. http://dept.stat.lsa.umich.edu/jerrick/courses/stat701/notes/mi.html.

  4. Graham, John W. 2009. Missing Data Analysis: Making It Work in the Real World. Annual Review of Psychology 60. Annual Reviews: 549-76.

  5. Rubin, Donald B. 1987. Multiple Imputation for Nonresponse in Surveys. Vol. 81. John Wiley & Sons.

*PS: Vielen Dank an Axel Klenk für den Tipp.

How to calculate Odds Ratios in Case-Control Studies using R

Intro

In June 2017 I've started working at the Clinical Trial Centre Leipzig at Leipzig University. Since my knowledge in statistics is rather poor, my employer offered me to attend some seminars in Medical Biometry at the University of Heidelberg. The first seminar I attended was called “Basics of Epidemiology”. At the first day, we learned how to calculate so called odds ratios in case-control studies using a simple pocket calculator.

In this blog post, I will show, how to calculate a simple odds ratio with 95% CI using R.

Data simulation

The data I'm using in this blog post were simulated using the wakefield package. The following code returns a data frame with 2 binary variables (Exposition and Disease) and 1.000 cases.

library(wakefield)

mydata <- data.frame(Exposition = group(n = 1000, x = c('yes', 'no'), 
                             prob = c(0.75, 0.25)),
                     Disease = group(n = 1000, x = c('yes', 'no'), 
                             prob = c(0.75, 0.25)))
dim(mydata)
## [1] 1000    2
head(mydata)
##   Exposition Disease
## 1        yes     yes
## 2         no      no
## 3         no     yes
## 4        yes     yes
## 5        yes     yes
## 6        yes     yes

Based on this data frame, we calculate a table showing how many patients with exposition vs. no exposition developed a disease vs. no disease.

tab <-table(mydata$Exposition, mydata$Disease)
tab
##      
##       yes  no
##   yes 569 210
##   no  163  58

Odds Ratio Calculation

In order to get to know whether the risk for developing a disease is significantly higher in patients having a certain exposition, we need to calculate the odds ratio and its 95% CI.

The following function will return a data frame containing these values.

# return odds ratio with 95%ci
f <- function(x) {
  or <- round((x[1] * x[4]) / (x[2] * x[3]), 2)
  cil <- round(exp(log(or) - 1.96 * sqrt(1/x[1] + 1/x[2] + 1/x[3] + 1/x[4])), 2)
  ciu <- round(exp(log(or) + 1.96 * sqrt(1/x[1] + 1/x[2] + 1/x[3] + 1/x[4])), 2)
  df <- data.frame(matrix(ncol = 3, nrow = 1, 
                              dimnames = list(NULL, c('CI_95_lower', 'OR', 'CI_95_upper'))))
  df[1,] <- rbind(c(cil, or, ciu))
  df <- as.data.frame(df)
}

Now, we can deploy the function on our table tab.

df.or <- f(tab)
knitr::kable(df.or, align = 'c')
CI_95_lower OR CI_95_upper
0.68 0.96 1.35

As the results indicate, patients with a disposition have no higher risk to develop a disease than patients having no disposition.

How to use R for matching samples (propensity score)

According to Wikipedia, propensity score matching (PSM) is a “statistical matching technique that attempts to estimate the effect of a treatment, policy, or other intervention by accounting for the covariates that predict receiving the treatment”. In a broader sense, propensity score analysis assumes that an unbiased comparison between samples can only be made when the subjects of both samples have similar characteristics. Thus, PSM can not only be used as “an alternative method to estimate the effect of receiving treatment when random assignment of treatments to subjects is not feasible” (Thavaneswaran 2008). It can also be used for the comparison of samples in epidemiological studies. Let's give an example:

Health-related quality of life (HRQOL) is considered an important outcome in cancer therapy. One of the most frequently used instruments to measure HRQOL in cancer patients is the core quality-of-life questionnaire of the European Organisation for Research and Treatment of Cancer. The EORTC QLQ-C30 is a 30-item instrument comprised of five functioning scales, nine symptom scales and one scale measuring Global quality of life. All scales have a score range between 0 and 100. While high scores of the symptom scales indicate a high burden of symptoms, high scores of the functioning scales and on the GQoL scale indicate better functioning resp. quality of life.

However, without having any reference point, it is difficult if not impossible to interpret the scores. Fortunately, the EORTC QLQ-C30 questionnaire was used in several general population surveys. Therefore, patient scores may be compared against scores of the general population. This makes it far easier to decide whether the burden of symptoms or functional impairments can be attributed to cancer (treatment) or not. PSM can be used to make both patient and population samples comparable by matching for relevant demographic characteristics like age and sex.

In this blog post, I show how to do PSM using R. A more comprehensive PSM guide can be found under: A Step-by-Step Guide to Propensity Score Matching in R.

Creating two random dataframes

Since we don't want to use real-world data in this blog post, we need to emulate the data. This can be easily done using the Wakefield package.

In a first step, we create a dataframe named df.patients. We want the dataframe to contain specifications of age and sex for 250 patients. The patients' age shall be between 30 and 78 years. Furthermore, 70% of patients shall be male.

set.seed(1234)
df.patients <- r_data_frame(n = 250, 
                            age(x = 30:78, 
                                name = 'Age'), 
                            sex(x = c("Male", "Female"), 
                                prob = c(0.70, 0.30), 
                                name = "Sex"))
df.patients$Sample <- as.factor('Patients')

The summary-function returns some basic information about the dataframe created. As we can see, the mean age of the patient sample is 53.7 and roughly 70% of the patients are male (69.2%).

summary(df.patients)
##       Age            Sex           Sample   
##  Min.   :30.00   Male  :173   Patients:250  
##  1st Qu.:42.00   Female: 77                 
##  Median :54.00                              
##  Mean   :53.71                              
##  3rd Qu.:66.00                              
##  Max.   :78.00

In a second step, we create another dataframe named df.population. We want this dataframe to comprise the same variables as df.patients with different specifications. With 18 to 80, the age-range of the population shall be wider than in the patient sample and the proportion of female and male patients shall be the same.

set.seed(1234)
df.population <- r_data_frame(n = 1000, 
                              age(x = 18:80, 
                                  name = 'Age'), 
                              sex(x = c("Male", "Female"), 
                                  prob = c(0.50, 0.50), 
                                  name = "Sex"))
df.population$Sample <- as.factor('Population')

The following table shows the sample's mean age (49.5 years) and the proportion of men (48.5%) and women (51.5%).

summary(df.population)
##       Age            Sex             Sample    
##  Min.   :18.00   Male  :485   Population:1000  
##  1st Qu.:34.00   Female:515                    
##  Median :50.00                                 
##  Mean   :49.46                                 
##  3rd Qu.:65.00                                 
##  Max.   :80.00

Merging the dataframes

Before we match the samples, we need to merge both dataframes. Based on the variable Sample, we create a new variable named Group (type logic) and a further variable (Distress) containing information about the individuals' level of distress. The Distress variable is created using the age-function of the Wakefield package. As we can see, women will have higher levels of distress.

mydata <- rbind(df.patients, df.population)
mydata$Group <- as.logical(mydata$Sample == 'Patients')
mydata$Distress <- ifelse(mydata$Sex == 'Male', age(nrow(mydata), x = 0:42, name = 'Distress'),
                                                age(nrow(mydata), x = 15:42, name = 'Distress'))

When we compare the distribution of age and sex in both samples, we discover significant differences:

pacman::p_load(tableone)
table1 <- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), 
                         data = mydata, 
                         factorVars = 'Sex', 
                         strata = 'Sample')
table1 <- print(table1, 
                printToggle = FALSE, 
                noSpaces = TRUE)
kable(table1[,1:3],  
      align = 'c', 
      caption = 'Table 1: Comparison of unmatched samples')
Patients Population p
n 250 1000
Age (mean (sd)) 53.71 (13.88) 49.46 (18.33) 0.001
Sex = Female (%) 77 (30.8) 515 (51.5) <0.001
Distress (mean (sd)) 22.86 (11.38) 25.13 (11.11) 0.004

Furthermore, the level of distress seems to be significantly higher in the population sample.

Matching the samples

Now, that we have completed preparation and inspection of data, we are going to match the two samples using the matchit-function of the MatchIt package. The method command method="nearest" specifies that the nearest neighbors method will be used. Other matching methods are exact matching, subclassification, optimal matching, genetic matching, and full matching (method = c("exact", "subclass", "optimal", ""genetic", "full")). The ratio command ratio = 1 indicates a one-to-one matching approach. With regard to our example, for each case in the patient sample exactly one case in the population sample will be matched. Please also note that the Group variable needs to be logic (TRUE vs. FALSE).

set.seed(1234)
match.it <- matchit(Group ~ Age + Sex, data = mydata, method="nearest", ratio=1)
a <- summary(match.it)

For further data presentation, we save the output of the summary-function into a variable named a.

After matching the samples, the size of the population sample was reduced to the size of the patient sample (n=250; see table 2).

kable(a$nn, digits = 2, align = 'c', 
      caption = 'Table 2: Sample sizes')
Control Treated
All 1000 250
Matched 250 250
Unmatched 750 0
Discarded 0 0

The following output shows, that the distributions of the variables Age and Sex are nearly identical after matching.

kable(a$sum.matched[c(1,2,4)], digits = 2, align = 'c', 
      caption = 'Table 3: Summary of balance for matched data')
Means Treated Means Control Mean Diff
distance 0.23 0.23 0.00
Age 53.71 53.65 0.06
SexMale 0.69 0.69 0.00
SexFemale 0.31 0.31 0.00

The distributions of propensity scores can be visualized using the plot-function which is part of the MatchIt package .

plot(match.it, type = 'jitter', interactive = FALSE)

plot of chunk plot

Saving the matched samples

Finally, the matched samples will be saved into a new dataframe named df.match.

df.match <- match.data(match.it)[1:ncol(mydata)]
rm(df.patients, df.population)

Eventually, we can check whether the differences in the level of distress between both samples are still significant.

pacman::p_load(tableone)
table4 <- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), 
                         data = df.match, 
                         factorVars = 'Sex', 
                         strata = 'Sample')
table4 <- print(table4, 
                printToggle = FALSE, 
                noSpaces = TRUE)
kable(table4[,1:3],  
      align = 'c', 
      caption = 'Table 4: Comparison of matched samples')
Patients Population p
n 250 250
Age (mean (sd)) 53.71 (13.88) 53.65 (13.86) 0.961
Sex = Female (%) 77 (30.8) 77 (30.8) 1.000
Distress (mean (sd)) 22.86 (11.38) 24.13 (11.88) 0.222

With a p-value of 0.222, Student's t-test does not indicate significant differences anymore. Thus, PSM helped to avoid an alpha mistake.


PS 1: The packages used in this blog post can be loaded/installed using the following code:

pacman::p_load(knitr, wakefield, MatchIt, tableone, captioner)

PS 2: Thanks very much to my colleague Katharina Kuba for for telling me about the MatchIt package.

Design a site like this with WordPress.com
Get started