New to R? Before you start, install R and RStudio from posit.co and work through the PREDICT Practical 1 (Getting started with R) and the Introduction to R: the chapters “Start”, “Set-up” (R, RStudio) and “The basics”, and the exercises at the bottom of the pages Creating data, Using data, Packages and External data.
In this tutorial you will go from raw data to a final prediction model for a time-to-event outcome. It has four parts:
Load data and make Table 1: describe the study population.
Describe the survival data: define the outcome, make a Kaplan-Meier curve, choose a time horizon and calculate the median follow-up time.
Multiple imputation: deal with missing data.
Fit the models: check non-linear effects, select predictors and present the final model.
Tutorial 2 continues with the final model and validates it
(shrinkage, optimism-corrected C-index, calibration and decision
curves). We use the pbc data (primary biliary cirrhosis)
from the survival package as dummy data.
Do not just read this page: run the code yourself.
You only learn R by doing it. Download the R Markdown file of this
tutorial (Tutorial_1_Develop_Prediction_Model.Rmd) by
clicking the Code button at the top right of this page
and choosing Download Rmd, and run every chunk in
RStudio, one by one. Read the output, change the code, break it and fix
it, and try the exercises before you open the answers. Reading the
results on this page is not a substitute.
To make sure everything works:
Create a folder for the course, for example
Predictive_models/, and put the downloaded
.Rmd file in it.
Open the .Rmd file in RStudio. When
you knit or run chunks from an Rmd file, R uses the folder that contains
the file as the working directory. Check this by running:
getwd()
The output should be the path of the folder from step 1. If it is
not, set it with setwd("path/to/your/folder") or in RStudio
via Session > Set Working Directory > To Source File
Location.
Do not use file paths in the code. All files are read from and saved to the working directory, so keep all files in this one folder.
Running the tutorial saves its output (for example the
Kaplan-Meier figure KM_plot_10_years.png and Table 1 as a
Word file) in this same folder.
The folder will then look like this:
Predictive_models/
├── Tutorial_1_Develop_Prediction_Model.Rmd
├── Tutorial_2_Validate_Prediction_Model.Rmd
├── Descriptives.docx (created by Tutorial 1)
└── KM_plot_10_years.png (created by Tutorial 1)
Start with a clean environment so that results do not depend on
objects left over from earlier work. A package is a collection
of functions that someone wrote and shared.
library(package) loads a package, so that you can use its
functions directly. We only load rms, because we use its
functions (such as cph, Predict and
fastbw) without the package:: prefix. All
other packages are called with the :: notation (see below),
so they do not have to be loaded, but they do have to be installed once.
We use:
survival: the data and survival analysisdplyr: data manipulationarsenal and flextable: Table 1 and saving
it to Wordsurvminer and prodlim: Kaplan-Meier curves
and follow-up timevisdat and mice: visualizing and imputing
missing datarms and Hmisc: fitting and presenting
prediction models (cph, fit.mult.impute)ggplot2 and ggpubr: saving and combining
figures# clear the environment: remove all objects that R currently remembers
rm(list=ls())
# load the rms package: after this, all its functions can be used directly
library(rms)
# The other packages must be installed once on your computer (not run here):
# install.packages(c("survival", "dplyr", "arsenal", "flextable", "survminer",
# "prodlim", "visdat", "mice", "Hmisc", "ggplot2", "ggpubr"))
:: notationIn this tutorial you will often see code like
dplyr::select(...). The notation
package::function() means “use the function
function from the package package”.
You do not have to write it this way once you have loaded a package
with library(): after library(dplyr),
select(...) works as well. We still write it out because it
is nice to see where each function comes from, and it avoids mistakes
when two packages have a function with the same name (R then uses the
one from the package that was loaded last). For a package that you do
not load with library(), such as openxlsx or
foreign in the templates below, the ::
notation lets you use one of its functions without loading the whole
package.
Some steps in this tutorial use random numbers, most importantly
multiple imputation. Without a seed, you would get different imputed
values, and therefore slightly different results, every time you run the
code. set.seed() fixes the starting point of R’s random
number generator, so everyone who runs the code obtains the same
numbers:
# the same seed gives the same random numbers
set.seed(1)
rnorm(3)
## [1] -0.6264538 0.1836433 -0.8356286
set.seed(1)
rnorm(3)
## [1] -0.6264538 0.1836433 -0.8356286
We set the random seed once, here at the top. In
addition, mice() has its own seed argument,
which we use when we impute, so that the imputations do not depend on
what was run before. If you change the order of the code or add extra
random steps, the numbers after that point can change.
# set the seed once for the whole tutorial
set.seed(1)
Everything in this tutorial is saved in one folder: the folder that contains this Rmd file. When you knit the document (or run the chunks in RStudio), the working directory is this folder, so tables, figures and R data files are saved right next to the Rmd and you do not need to set any file paths.
Your turn: check where your files will be saved.
# all files are saved in the working directory, which is the folder of this Rmd
getwd()
Every prediction model project starts by looking at the data. In this part you load the data, decide which variables are candidate predictors, make sure every variable has the correct type, and create a “Table 1” that describes your study population. You save this table to Word so you can use it in a manuscript.
In your own project you will read your own data into R, most often
from an Excel or SPSS file. The next two blocks show how to do that.
They are not run in this tutorial
(eval=FALSE), because we do not have your files; they are
here as templates for your own projects. After them, we load the dummy
data that we actually use.
openxlsx::read.xlsx() reads an .xlsx file
into a data frame (older .xls files must first be saved as
.xlsx). Replace the path with the location of your file. If
your file is in the same folder as this Rmd, the file name alone is
enough. If the data are not on the first sheet, add the argument
sheet="Sheet name" (or the sheet number). Install the
package once with install.packages("openxlsx") if you do
not have it yet.
# load Excel data (this block is not run)
data <- openxlsx::read.xlsx("C:/.../file_data_name.xlsx")
foreign::read.spss() reads an SPSS .sav
file. With to.data.frame=TRUE the result is a data frame
instead of a list. SPSS stores dates as seconds since a starting date,
so origin="01-01-1970" tells R how to convert them to
dates. Variables with value labels in SPSS become factors in R, so check
afterwards that every variable has the type you expect (see the next
section).
# load SPSS file (this block is not run)
data <- foreign::read.spss(file="C:/.../file_data_name.sav", to.data.frame=TRUE, origin="01-01-1970")
Here we use the dummy data from the survival package.
The data dictionary is available at https://cran.r-project.org/web/packages/survival/survival.pdf.
# load the pbc data set from the survival package into R's memory
data(pbc, package="survival")
# copy it to an object with the general name `data` (<- assigns a value to a name)
data <- pbc
We also add a center variable. Note: this
variable is artificial. It simply alternates between 1 and 2,
so it carries no real between-center differences. That is fine for a
demonstration; we use it in Tutorial 2 to illustrate internal-external
cross-validation.
# add center variable for illustration (artificial, fine for a demo)
# rep(1:2, length.out=n) repeats 1, 2, 1, 2, ... until we have n values, with
# n = nrow(data) the number of patients; `data$center <-` stores them as a new column
data$center <- rep(1:2, length.out=nrow(data))
Picking a column: $ versus
[[. Here and below we use data$center
to pick a column, which is short to type. R has two ways to take one
column out of a data frame, and they behave differently:
[ keeps the container, [[ takes one
element out of it: data["center"] is still a data frame,
whereas data[["center"]] is the column itself.
data$center also gives the column itself.$ does partial matching: data$cen
silently gives you the center column. In the console that
is convenient, but in a script a typo returns the wrong column instead
of an error.$ takes the name literally, so it cannot use a column
name that is stored in a variable. [[ evaluates the
variable first, which is what makes loops over variable names possible
(we use this in Part 3).# the column itself
head(data$center)
## [1] 1 2 1 2 1 2
head(data[["center"]])
## [1] 1 2 1 2 1 2
class(data["center"]) # a data frame with one column
## [1] "data.frame"
class(data[["center"]]) # the column itself
## [1] "integer"
# partial matching with $
head(data$cen)
## [1] 1 2 1 2 1 2
# column name stored in a variable
my_var <- "center"
data$my_var # NULL: looks for a column called my_var
## NULL
head(data[[my_var]]) # works
## [1] 1 2 1 2 1 2
In scripts, [["name"]] is therefore safer, and it is the
only option when the column name is stored in a variable.
Select the outcome variables (time and
status) and the variables that you want to describe and use
as candidate predictors, based on literature and clinical knowledge.
Binary and categorical variables must be factor variables,
otherwise R treats them as numbers and reports a mean and SD for
them.
Two variables deserve attention: trt (treatment, coded
1/2) and edema (coded 0, 0.5, 1) are categories, not
quantities, so they are also converted to factors. We set every factor
separately with factor(): levels are the
values as they are stored in the data and labels are the
names that we want to see in tables, figures and model output. The
labels must be in the same order as the levels. Check the codebook of
your data to see what each value means (in pbc,
edema 0 = no edema, 0.5 = edema untreated or successfully
treated, 1 = edema despite diuretic therapy, and sex is
stored as “m” and “f”).
Note: check how binary variables are coded.
trt is coded 1/2 here. This often happens when you load
data from SPSS, where categories are stored as 1/2 instead of 0/1.
Always check this (for example with table(data$trt) or
str(data)), because a variable coded 1/2 is easily mistaken
for a quantity, and the reference category is not obvious. It is
probably best to recode binary variables to 0/1 (for example 0 =
no/control, 1 = yes/treatment), or to give them readable labels as we do
below. Never rely on as.factor() alone, because it keeps
the stored values (1 and 2) as labels.
# select the outcome and candidate variables based on literature
# (dplyr::select() keeps only the columns that we name)
data <- dplyr::select(data, time, status,
trt, age, sex, ascites, hepato, spiders, edema,
bili, chol, albumin, copper, alk.phos, ast, trig,
platelet, protime, stage, center)
# binary and categorical variables should be factor variables:
# factor(x, levels=..., labels=...) takes the values that are stored in the data
# (levels) and gives each of them a readable label (in the same order)
data$trt <- factor(data$trt, levels=c(1, 2),
labels=c("D-penicillamine", "Placebo"))
data$sex <- factor(data$sex, levels=c("m", "f"),
labels=c("Male", "Female"))
data$ascites <- factor(data$ascites, levels=c(0, 1),
labels=c("No", "Yes"))
data$hepato <- factor(data$hepato, levels=c(0, 1),
labels=c("No", "Yes"))
data$spiders <- factor(data$spiders, levels=c(0, 1),
labels=c("No", "Yes"))
data$edema <- factor(data$edema, levels=c(0, 0.5, 1),
labels=c("No edema", "Untreated or successfully treated",
"Despite diuretics"))
data$stage <- factor(data$stage, levels=1:4,
labels=c("1", "2", "3", "4"))
data$center <- factor(data$center, levels=c(1, 2),
labels=c("Center 1", "Center 2"))
Your turn: run str(data) and check that
every variable has the type you expect.
Check the labels. factor() silently
turns every value that is not listed in levels into a
missing value (NA), for example when a level is misspelled
or the data are coded differently than you think. After setting the
labels, always count the values per level and include the missing values
with useNA="ifany", and compare the counts with those
before recoding:
# table() counts how often each value occurs; useNA="ifany" also counts the
# missing values (NA), so that wrongly labelled values show up
table(data$trt, useNA="ifany")
##
## D-penicillamine Placebo <NA>
## 158 154 106
table(data$sex, useNA="ifany")
##
## Male Female
## 44 374
table(data$edema, useNA="ifany")
##
## No edema Untreated or successfully treated
## 354 44
## Despite diuretics
## 20
The first level is the reference category in models:
the coefficients of a factor compare each other level with the first
one. That is why the model output later shows names such as
ascites=Yes (compared with “No”) and
sex=Female (compared with “Male”). Use
relevel(x, ref="...") to choose another reference
category.
arsenal::tableby() summarizes each variable, here by
center. We show the mean (SD), median (Q1, Q3), range and
number of missing values for continuous variables, and counts
(percentages) and missing values for categorical variables. The outcome
variables time and status are left out of the
table. We switch off the statistical tests (test=FALSE)
because a Table 1 in a prediction model paper describes the population
and does not compare groups.
# center ~ . means: one column per center, describing all other variables (.)
# data=...: all variables except time and status (-time, -status leaves them out)
table_one <- arsenal::tableby(center ~ .,
data=dplyr::select(data, -time, -status),
total.pos="before", # put the total column first
digits=1, # number of decimals
test=FALSE, # no p-values
total=TRUE, # add a column for all patients
numeric.test="kwt",
cat.test="chisq",
cat.simplify=TRUE, # one row for binary variables
# statistics shown for continuous variables ...
numeric.stats=c("meansd", "medianq1q3",
"range", "Nmiss2"),
# ... and for categorical variables
cat.stats=c("countpct", "Nmiss2"),
# readable names for the statistics
stats.labels=list(
meansd="Mean (SD)",
medianq1q3="Median (Q1, Q3)",
range="Min - Max",
Nmiss2="Missing"))
Here is the table as it will look in the report:
summary(table_one, text="html")
| Total (N=418) | Center 1 (N=209) | Center 2 (N=209) | |
|---|---|---|---|
| trt | |||
| D-penicillamine | 158 (50.6%) | 71 (45.5%) | 87 (55.8%) |
| Placebo | 154 (49.4%) | 85 (54.5%) | 69 (44.2%) |
| Missing | 106 | 53 | 53 |
| age | |||
| Mean (SD) | 50.7 (10.4) | 50.5 (10.0) | 51.0 (10.9) |
| Median (Q1, Q3) | 51.0 (42.8, 58.2) | 49.6 (43.0, 58.0) | 52.0 (42.8, 58.4) |
| Min - Max | 26.3 - 78.4 | 29.6 - 78.4 | 26.3 - 76.7 |
| Missing | 0 | 0 | 0 |
| sex | |||
| Male | 44 (10.5%) | 20 (9.6%) | 24 (11.5%) |
| Female | 374 (89.5%) | 189 (90.4%) | 185 (88.5%) |
| Missing | 0 | 0 | 0 |
| ascites | |||
| No | 288 (92.3%) | 141 (90.4%) | 147 (94.2%) |
| Yes | 24 (7.7%) | 15 (9.6%) | 9 (5.8%) |
| Missing | 106 | 53 | 53 |
| hepato | |||
| No | 152 (48.7%) | 76 (48.7%) | 76 (48.7%) |
| Yes | 160 (51.3%) | 80 (51.3%) | 80 (51.3%) |
| Missing | 106 | 53 | 53 |
| spiders | |||
| No | 222 (71.2%) | 111 (71.2%) | 111 (71.2%) |
| Yes | 90 (28.8%) | 45 (28.8%) | 45 (28.8%) |
| Missing | 106 | 53 | 53 |
| edema | |||
| No edema | 354 (84.7%) | 172 (82.3%) | 182 (87.1%) |
| Untreated or successfully treated | 44 (10.5%) | 25 (12.0%) | 19 (9.1%) |
| Despite diuretics | 20 (4.8%) | 12 (5.7%) | 8 (3.8%) |
| Missing | 0 | 0 | 0 |
| bili | |||
| Mean (SD) | 3.2 (4.4) | 3.2 (4.5) | 3.2 (4.3) |
| Median (Q1, Q3) | 1.4 (0.8, 3.4) | 1.3 (0.8, 3.2) | 1.5 (0.8, 3.5) |
| Min - Max | 0.3 - 28.0 | 0.3 - 24.5 | 0.3 - 28.0 |
| Missing | 0 | 0 | 0 |
| chol | |||
| Mean (SD) | 369.5 (231.9) | 368.9 (221.7) | 370.2 (243.2) |
| Median (Q1, Q3) | 309.5 (249.5, 400.0) | 316.0 (254.5, 400.0) | 302.0 (246.0, 396.0) |
| Min - Max | 120.0 - 1775.0 | 120.0 - 1712.0 | 127.0 - 1775.0 |
| Missing | 134 | 62 | 72 |
| albumin | |||
| Mean (SD) | 3.5 (0.4) | 3.5 (0.4) | 3.5 (0.4) |
| Median (Q1, Q3) | 3.5 (3.2, 3.8) | 3.6 (3.2, 3.8) | 3.5 (3.3, 3.8) |
| Min - Max | 2.0 - 4.6 | 2.0 - 4.5 | 2.2 - 4.6 |
| Missing | 0 | 0 | 0 |
| copper | |||
| Mean (SD) | 97.6 (85.6) | 99.0 (81.9) | 96.3 (89.5) |
| Median (Q1, Q3) | 73.0 (41.2, 123.0) | 75.0 (44.0, 129.0) | 67.0 (39.0, 121.0) |
| Min - Max | 4.0 - 588.0 | 4.0 - 558.0 | 11.0 - 588.0 |
| Missing | 108 | 53 | 55 |
| alk.phos | |||
| Mean (SD) | 1982.7 (2140.4) | 2022.6 (2118.2) | 1942.8 (2168.5) |
| Median (Q1, Q3) | 1259.0 (871.5, 1980.0) | 1300.0 (857.8, 2040.5) | 1224.0 (877.5, 1929.2) |
| Min - Max | 289.0 - 13862.4 | 377.0 - 12258.8 | 289.0 - 13862.4 |
| Missing | 106 | 53 | 53 |
| ast | |||
| Mean (SD) | 122.6 (56.7) | 123.3 (56.3) | 121.9 (57.3) |
| Median (Q1, Q3) | 114.7 (80.6, 151.9) | 110.5 (82.2, 153.4) | 116.6 (80.6, 148.0) |
| Min - Max | 26.4 - 457.2 | 41.9 - 338.0 | 26.4 - 457.2 |
| Missing | 106 | 53 | 53 |
| trig | |||
| Mean (SD) | 124.7 (65.1) | 125.4 (68.6) | 123.9 (61.5) |
| Median (Q1, Q3) | 108.0 (84.2, 151.0) | 112.0 (85.0, 149.0) | 106.0 (84.0, 153.0) |
| Min - Max | 33.0 - 598.0 | 44.0 - 598.0 | 33.0 - 382.0 |
| Missing | 136 | 64 | 72 |
| platelet | |||
| Mean (SD) | 257.0 (98.3) | 260.2 (98.6) | 253.8 (98.2) |
| Median (Q1, Q3) | 251.0 (188.5, 318.0) | 256.0 (195.0, 319.0) | 244.0 (181.0, 317.2) |
| Min - Max | 62.0 - 721.0 | 62.0 - 563.0 | 70.0 - 721.0 |
| Missing | 11 | 6 | 5 |
| protime | |||
| Mean (SD) | 10.7 (1.0) | 10.8 (1.1) | 10.7 (0.9) |
| Median (Q1, Q3) | 10.6 (10.0, 11.1) | 10.6 (10.1, 11.1) | 10.6 (10.0, 11.0) |
| Min - Max | 9.0 - 18.0 | 9.3 - 18.0 | 9.0 - 14.1 |
| Missing | 2 | 1 | 1 |
| stage | |||
| 1 | 21 (5.1%) | 9 (4.4%) | 12 (5.8%) |
| 2 | 92 (22.3%) | 47 (22.9%) | 45 (21.7%) |
| 3 | 155 (37.6%) | 67 (32.7%) | 88 (42.5%) |
| 4 | 144 (35.0%) | 82 (40.0%) | 62 (30.0%) |
| Missing | 6 | 4 | 2 |
# convert the table to a data frame, make a flextable of it and save it as a
# Word file in the working directory. The pipe |> passes the result on to the
# next function
flextable::flextable(data.frame(summary(table_one, text=TRUE))) |>
flextable::save_as_docx(path="Descriptives.docx")
Try each exercise yourself first. Then open the field below it to see the answer, the code and its results. The answers use new object names, so they do not change anything in the tutorial.
Exercise 1. Which variables have the most missing values? Think about what this means for the modelling in Part 3.
Count the missing values per variable, and the percentage of patients with a complete record:
# number of missing values per variable, most missing first
sort(colSums(is.na(data)), decreasing=TRUE)
## trig chol copper trt ascites hepato spiders alk.phos
## 136 134 108 106 106 106 106 106
## ast platelet stage protime time status age sex
## 106 11 6 2 0 0 0 0
## edema bili albumin center
## 0 0 0 0
# percentage of patients without any missing value
round(100*mean(complete.cases(data)), 1)
## [1] 66
Several variables are missing for many patients. If you deleted every
patient with at least one missing value (complete case analysis), you
would throw away a large part of the data. In pbc, the
patients who did not take part in the randomized trial (missing
trt) also have many other missing values, so the
missingness is related to other variables. Both are reasons to impute
(Part 3) instead of using complete cases.
Exercise 2. Change center ~ . to
trt ~ . or sex ~ . and look at the new table.
Which would be most useful in a paper?
table_trt <- arsenal::tableby(trt ~ .,
data=dplyr::select(data, -time, -status),
test=FALSE, digits=1)
summary(table_trt, text="html")
| D-penicillamine (N=158) | Placebo (N=154) | Total (N=312) | |
|---|---|---|---|
| age | |||
| Mean (SD) | 51.4 (11.0) | 48.6 (10.0) | 50.0 (10.6) |
| Range | 26.3 - 78.4 | 30.6 - 74.5 | 26.3 - 78.4 |
| sex | |||
| Male | 21 (13.3%) | 15 (9.7%) | 36 (11.5%) |
| Female | 137 (86.7%) | 139 (90.3%) | 276 (88.5%) |
| ascites | |||
| No | 144 (91.1%) | 144 (93.5%) | 288 (92.3%) |
| Yes | 14 (8.9%) | 10 (6.5%) | 24 (7.7%) |
| hepato | |||
| No | 85 (53.8%) | 67 (43.5%) | 152 (48.7%) |
| Yes | 73 (46.2%) | 87 (56.5%) | 160 (51.3%) |
| spiders | |||
| No | 113 (71.5%) | 109 (70.8%) | 222 (71.2%) |
| Yes | 45 (28.5%) | 45 (29.2%) | 90 (28.8%) |
| edema | |||
| No edema | 132 (83.5%) | 131 (85.1%) | 263 (84.3%) |
| Untreated or successfully treated | 16 (10.1%) | 13 (8.4%) | 29 (9.3%) |
| Despite diuretics | 10 (6.3%) | 10 (6.5%) | 20 (6.4%) |
| bili | |||
| Mean (SD) | 2.9 (3.6) | 3.6 (5.3) | 3.3 (4.5) |
| Range | 0.3 - 20.0 | 0.3 - 28.0 | 0.3 - 28.0 |
| chol | |||
| N-Miss | 18 | 10 | 28 |
| Mean (SD) | 365.0 (209.5) | 373.9 (252.5) | 369.5 (231.9) |
| Range | 127.0 - 1712.0 | 120.0 - 1775.0 | 120.0 - 1775.0 |
| albumin | |||
| Mean (SD) | 3.5 (0.4) | 3.5 (0.4) | 3.5 (0.4) |
| Range | 2.1 - 4.6 | 2.0 - 4.4 | 2.0 - 4.6 |
| copper | |||
| N-Miss | 1 | 1 | 2 |
| Mean (SD) | 97.6 (90.6) | 97.7 (80.5) | 97.6 (85.6) |
| Range | 9.0 - 588.0 | 4.0 - 558.0 | 4.0 - 588.0 |
| alk.phos | |||
| Mean (SD) | 2021.3 (2183.4) | 1943.0 (2101.7) | 1982.7 (2140.4) |
| Range | 369.0 - 11552.0 | 289.0 - 13862.4 | 289.0 - 13862.4 |
| ast | |||
| Mean (SD) | 120.2 (54.5) | 125.0 (58.9) | 122.6 (56.7) |
| Range | 26.4 - 338.0 | 28.4 - 457.2 | 26.4 - 457.2 |
| trig | |||
| N-Miss | 19 | 11 | 30 |
| Mean (SD) | 124.1 (71.5) | 125.3 (58.5) | 124.7 (65.1) |
| Range | 33.0 - 598.0 | 44.0 - 432.0 | 33.0 - 598.0 |
| platelet | |||
| N-Miss | 2 | 2 | 4 |
| Mean (SD) | 258.8 (100.3) | 265.2 (90.7) | 261.9 (95.6) |
| Range | 62.0 - 563.0 | 71.0 - 487.0 | 62.0 - 563.0 |
| protime | |||
| Mean (SD) | 10.7 (0.9) | 10.8 (1.1) | 10.7 (1.0) |
| Range | 9.0 - 14.1 | 9.2 - 17.1 | 9.0 - 17.1 |
| stage | |||
| 1 | 12 (7.6%) | 4 (2.6%) | 16 (5.1%) |
| 2 | 35 (22.2%) | 32 (20.8%) | 67 (21.5%) |
| 3 | 56 (35.4%) | 64 (41.6%) | 120 (38.5%) |
| 4 | 55 (34.8%) | 54 (35.1%) | 109 (34.9%) |
| center | |||
| Center 1 | 71 (44.9%) | 85 (55.2%) | 156 (50.0%) |
| Center 2 | 87 (55.1%) | 69 (44.8%) | 156 (50.0%) |
Which split is most useful depends on the purpose. In a prediction
model paper, Table 1 usually describes the whole development cohort,
sometimes split by center or by the outcome. A table by treatment helps
to check whether the treatment groups are comparable. Note that patients
with missing trt do not appear in either column.
Exercise 3. Remove "Nmiss2" from
cat.stats and see what changes.
table_no_miss <- arsenal::tableby(center ~ .,
data=dplyr::select(data, -time, -status),
total.pos="before",
digits=1,
test=FALSE,
total=TRUE,
cat.simplify=TRUE,
numeric.stats=c("meansd", "medianq1q3",
"range", "Nmiss2"),
cat.stats=c("countpct"),
stats.labels=list(
meansd="Mean (SD)",
medianq1q3="Median (Q1, Q3)",
range="Min - Max",
Nmiss2="Missing"))
summary(table_no_miss, text="html")
| Total (N=418) | Center 1 (N=209) | Center 2 (N=209) | |
|---|---|---|---|
| trt | 154 (49.4%) | 85 (54.5%) | 69 (44.2%) |
| age | |||
| Mean (SD) | 50.7 (10.4) | 50.5 (10.0) | 51.0 (10.9) |
| Median (Q1, Q3) | 51.0 (42.8, 58.2) | 49.6 (43.0, 58.0) | 52.0 (42.8, 58.4) |
| Min - Max | 26.3 - 78.4 | 29.6 - 78.4 | 26.3 - 76.7 |
| Missing | 0 | 0 | 0 |
| sex | 374 (89.5%) | 189 (90.4%) | 185 (88.5%) |
| ascites | 24 (7.7%) | 15 (9.6%) | 9 (5.8%) |
| hepato | 160 (51.3%) | 80 (51.3%) | 80 (51.3%) |
| spiders | 90 (28.8%) | 45 (28.8%) | 45 (28.8%) |
| edema | |||
| No edema | 354 (84.7%) | 172 (82.3%) | 182 (87.1%) |
| Untreated or successfully treated | 44 (10.5%) | 25 (12.0%) | 19 (9.1%) |
| Despite diuretics | 20 (4.8%) | 12 (5.7%) | 8 (3.8%) |
| bili | |||
| Mean (SD) | 3.2 (4.4) | 3.2 (4.5) | 3.2 (4.3) |
| Median (Q1, Q3) | 1.4 (0.8, 3.4) | 1.3 (0.8, 3.2) | 1.5 (0.8, 3.5) |
| Min - Max | 0.3 - 28.0 | 0.3 - 24.5 | 0.3 - 28.0 |
| Missing | 0 | 0 | 0 |
| chol | |||
| Mean (SD) | 369.5 (231.9) | 368.9 (221.7) | 370.2 (243.2) |
| Median (Q1, Q3) | 309.5 (249.5, 400.0) | 316.0 (254.5, 400.0) | 302.0 (246.0, 396.0) |
| Min - Max | 120.0 - 1775.0 | 120.0 - 1712.0 | 127.0 - 1775.0 |
| Missing | 134 | 62 | 72 |
| albumin | |||
| Mean (SD) | 3.5 (0.4) | 3.5 (0.4) | 3.5 (0.4) |
| Median (Q1, Q3) | 3.5 (3.2, 3.8) | 3.6 (3.2, 3.8) | 3.5 (3.3, 3.8) |
| Min - Max | 2.0 - 4.6 | 2.0 - 4.5 | 2.2 - 4.6 |
| Missing | 0 | 0 | 0 |
| copper | |||
| Mean (SD) | 97.6 (85.6) | 99.0 (81.9) | 96.3 (89.5) |
| Median (Q1, Q3) | 73.0 (41.2, 123.0) | 75.0 (44.0, 129.0) | 67.0 (39.0, 121.0) |
| Min - Max | 4.0 - 588.0 | 4.0 - 558.0 | 11.0 - 588.0 |
| Missing | 108 | 53 | 55 |
| alk.phos | |||
| Mean (SD) | 1982.7 (2140.4) | 2022.6 (2118.2) | 1942.8 (2168.5) |
| Median (Q1, Q3) | 1259.0 (871.5, 1980.0) | 1300.0 (857.8, 2040.5) | 1224.0 (877.5, 1929.2) |
| Min - Max | 289.0 - 13862.4 | 377.0 - 12258.8 | 289.0 - 13862.4 |
| Missing | 106 | 53 | 53 |
| ast | |||
| Mean (SD) | 122.6 (56.7) | 123.3 (56.3) | 121.9 (57.3) |
| Median (Q1, Q3) | 114.7 (80.6, 151.9) | 110.5 (82.2, 153.4) | 116.6 (80.6, 148.0) |
| Min - Max | 26.4 - 457.2 | 41.9 - 338.0 | 26.4 - 457.2 |
| Missing | 106 | 53 | 53 |
| trig | |||
| Mean (SD) | 124.7 (65.1) | 125.4 (68.6) | 123.9 (61.5) |
| Median (Q1, Q3) | 108.0 (84.2, 151.0) | 112.0 (85.0, 149.0) | 106.0 (84.0, 153.0) |
| Min - Max | 33.0 - 598.0 | 44.0 - 598.0 | 33.0 - 382.0 |
| Missing | 136 | 64 | 72 |
| platelet | |||
| Mean (SD) | 257.0 (98.3) | 260.2 (98.6) | 253.8 (98.2) |
| Median (Q1, Q3) | 251.0 (188.5, 318.0) | 256.0 (195.0, 319.0) | 244.0 (181.0, 317.2) |
| Min - Max | 62.0 - 721.0 | 62.0 - 563.0 | 70.0 - 721.0 |
| Missing | 11 | 6 | 5 |
| protime | |||
| Mean (SD) | 10.7 (1.0) | 10.8 (1.1) | 10.7 (0.9) |
| Median (Q1, Q3) | 10.6 (10.0, 11.1) | 10.6 (10.1, 11.1) | 10.6 (10.0, 11.0) |
| Min - Max | 9.0 - 18.0 | 9.3 - 18.0 | 9.0 - 14.1 |
| Missing | 2 | 1 | 1 |
| stage | |||
| 1 | 21 (5.1%) | 9 (4.4%) | 12 (5.8%) |
| 2 | 92 (22.3%) | 47 (22.9%) | 45 (21.7%) |
| 3 | 155 (37.6%) | 67 (32.7%) | 88 (42.5%) |
| 4 | 144 (35.0%) | 82 (40.0%) | 62 (30.0%) |
The “Missing” rows of the categorical variables disappear, so you can no longer see how many values are missing for them. Always report missing data in Table 1.
Exercise 4. Open Descriptives.docx in
your folder. Note that you still need to adjust things before you can
put this table in a manuscript. What would you change?
The Word table is created from the same table as above, with generic column names:
names(data.frame(summary(table_one, text=TRUE)))
## [1] "Var.1" "Total..N.418." "Center.1..N.209." "Center.2..N.209."
Things to check or improve in Word: the column headers (replace the generic names by clear ones), the number of decimals, readable variable labels and units (for example “alk.phos” is not a label for a paper), and add a table caption and a footnote that explains the abbreviations.
Our outcome is time to death, which is time-to-event data: some patients are followed until they die, others are censored (still alive at last contact). In this part you define the outcome, draw a Kaplan-Meier curve, restrict the outcome to a clinically relevant time horizon, and calculate the median follow-up time correctly.
Time-to-event data are the most common type of outcome data in
clinical prediction research, so it is worth learning how to work with
them. (If your own outcome is binary without follow-up time, you can
skip the parts about follow-up.) We create a Surv object
that combines follow-up time and event status. We convert time from days
to years, and dichotomize the event: in pbc the status is
coded 0 = censored, 1 = transplant, 2 = dead, and we model death.
Important: with status==2 we treat a
liver transplant as a censored observation. We do
not treat transplant as a competing risk. This is a
simplification: censoring at transplant assumes that patients who
receive a transplant would otherwise have the same risk of death as
those who remain in follow-up, which is unlikely to be true. If
transplant is an important competing event in your own project, use a
competing risk approach (for example cause-specific hazards or
Fine-Gray).
# generate the Surv object for the outcome: it combines follow-up time and event
# time: transform the time from days to years (divide by 365)
# event: status==2 is TRUE for patients who died and FALSE for the others;
# as.numeric() turns TRUE/FALSE into 1/0
# (transplant is treated as censored, NOT as a competing risk)
S <- survival::Surv(time=data$time/365, event=as.numeric(data$status==2))
The Kaplan-Meier method estimates survival over time while correctly
handling censoring. We make one curve per treatment group, here for all
follow-up (unrestricted, S).
# make Kaplan-Meier curve (all follow-up)
# the formula S ~ trt means: one survival curve for each treatment group
KM_unrestricted <- survival::survfit(S ~ trt, data=data)
Note: some patients are silently dropped.
trt is missing for some patients (those not in the
randomized trial), and survfit() removes them without a
warning. Check how many:
# number of patients with missing treatment, these are not in the curves
sum(is.na(data$trt))
## [1] 106
Only the remaining patients are used in the curves. In the prediction models in Part 4 we do keep everyone, because missing values are imputed.
First look at the curve with the basic plot()
function:
# basic Kaplan-Meier curve: one color per treatment group, with axis labels
plot(KM_unrestricted, col=c("blue", "red"), xlab="Time (years)", ylab="Survival")
For a nicer figure with confidence intervals and a risk table we use
ggsurvplot() from survminer.
# Kaplan-Meier curve with risk table (all follow-up)
KM_plot_unrestricted <- survminer::ggsurvplot(KM_unrestricted,
data=data,
risk.table=TRUE, # table below the curve
risk.table.y.text=FALSE, # no group names in the table
legend="top",
legend.title="Treatment",
legend.labs=levels(data$trt), # names of the groups
xlab="Time (years)",
conf.int=TRUE) # confidence intervals
# typing the name of an object prints it
KM_plot_unrestricted
Two details in the code make the figure tidier.
legend.labs=levels(data$trt) replaces the default labels
(trt=D-penicillamine, trt=Placebo) of the two
treatment arms in the legend by the factor labels we set earlier; to use
other names, give your own vector, in the same order as the levels, for
example legend.labs=c("Penicillamine", "Control"). And
risk.table.y.text=FALSE removes the long group names from
the y-axis of the risk table (the colors in the table match the legend),
so that the y-axes of the curve and the table line up.
Look at the right end of the curves: the confidence intervals become wider, and the risk table shows why: fewer and fewer patients are still followed after about 10 years, so the estimates there are less certain. A prediction model should not rely much on this part of the follow-up. Moreover, a prediction horizon of 10 years may be more clinically useful than one of 12 years. That is why we truncate the outcome at 10 years in the next section; for the purpose of this exercise we use 10 years as the horizon.
A prediction model predicts risk at a specific time point. We choose 10 years. Everyone followed longer than 10 years is censored at 10 years, because events after the horizon are not part of the outcome we predict.
# limit the outcome to your time horizon
horizon <- 10
S_10 <- S # start from a copy of the full outcome
# S_10[, 1] is the time column and S_10[, 2] is the event column;
# S_10[, 1] > horizon is TRUE for everyone followed longer than 10 years
S_10[S_10[, 1]>horizon, 2] <- 0 # for them: no event (censored)
S_10[S_10[, 1]>horizon, 1] <- horizon # and the time becomes 10 years
Now fit the Kaplan-Meier curve for the restricted outcome
S_10 and make the same figure as before:
# the same Kaplan-Meier fit, now for the outcome restricted to 10 years
KM_10 <- survival::survfit(S_10 ~ trt, data=data)
# Kaplan-Meier curve with risk table, restricted to 10 years
KM_plot <- survminer::ggsurvplot(KM_10,
data=data,
risk.table=TRUE,
risk.table.y.text=FALSE,
legend="top",
legend.title="Treatment",
legend.labs=levels(data$trt),
xlab="Time (years)",
conf.int=TRUE)
KM_plot
# combine the curve and the risk table and save
ggplot2::ggsave(filename="KM_plot_10_years.png",
plot=ggpubr::ggarrange(KM_plot$plot, KM_plot$table,
nrow=2, heights=c(3, 1), align="v"),
width=7, height=7, dpi=300)
Now display the Kaplan-Meier estimates at the horizon, for the unrestricted and the restricted outcome:
# display the Kaplan-Meier estimates (survival probability) at the horizon,
# time=horizon asks for the estimate at 10 years only
summary(KM_unrestricted, time=horizon)
## Call: survfit(formula = S ~ trt, data = data)
##
## 106 observations deleted due to missingness
## trt=D-penicillamine
## time n.risk n.event survival std.err lower 95% CI
## 10.0000 16.0000 63.0000 0.4247 0.0603 0.3216
## upper 95% CI
## 0.5611
##
## trt=Placebo
## time n.risk n.event survival std.err lower 95% CI
## 10.0000 16.0000 57.0000 0.4575 0.0612 0.3520
## upper 95% CI
## 0.5946
summary(KM_10, time=horizon)
## Call: survfit(formula = S_10 ~ trt, data = data)
##
## 106 observations deleted due to missingness
## trt=D-penicillamine
## time n.risk n.event survival std.err lower 95% CI
## 10.0000 16.0000 63.0000 0.4247 0.0603 0.3216
## upper 95% CI
## 0.5611
##
## trt=Placebo
## time n.risk n.event survival std.err lower 95% CI
## 10.0000 16.0000 57.0000 0.4575 0.0612 0.3520
## upper 95% CI
## 0.5946
The median of the observed times is not the median follow-up
time, because patients who died have a short observed time and drag the
median down. The correct estimate uses the reverse Kaplan-Meier
method: we flip the roles of event and censoring, so the “event” is
being censored and we estimate how long patients are followed in a
typical situation. In survfit() this means that the event
indicator becomes 1 - status.
First plot the reverse Kaplan-Meier curve. It shows the probability that a patient is still being followed at each time point; the time at which this curve crosses 0.5 is the median follow-up time.
# reverse Kaplan-Meier: the "event" is being censored (not dying)
# S[, "time"] is the time column of S; 1 - S[, "status"] flips 1 and 0, so
# that being censored counts as the event
reverse_KM <- survival::survfit(survival::Surv(S[, "time"], 1 - S[, "status"]) ~ 1,
data=data, conf.type="plain")
survminer::ggsurvplot(reverse_KM,
data=data,
conf.int=TRUE,
legend="none",
xlab="Time (years)",
ylab="Probability of still being followed")
Then calculate the median follow-up time from this curve:
# median follow-up time with the reverse Kaplan-Meier method (in years)
# probs=0.5 asks for the 50th percentile, the median
stats::quantile(reverse_KM, probs=0.5)
## $quantile
## 50
## 6.479452
##
## $lower
## 50
## 6.093151
##
## $upper
## 50
## 7.021918
The result gives the median follow-up time (with
probs=0.5) with its confidence interval. We use
conf.type="plain" in survfit() to get the same
confidence interval as prodlim below; the default of
survfit() is conf.type="log", which gives a
slightly different lower limit.
The prodlim package does the same in one line: with
reverse=TRUE it flips event and censoring for you, so you
do not have to define 1 - status yourself. Install it once
with install.packages("prodlim") if you do not have it yet.
The row for 0.5 gives the same median follow-up time and
confidence interval. Do not confuse this interval with the line “Median
with interquartile range (IQR)” below the table: that is the median
together with the 25th and 75th percentile of the follow-up times (the
rows for 0.75 and 0.25), which describes how
much the follow-up varies between patients. It is not a confidence
interval of the median.
# the same in one line with prodlim
stats::quantile(prodlim::prodlim(prodlim::Hist(time/365, status==2) ~ 1,
data=data, reverse=TRUE))
## Quantiles of the potential follow up time distribution based on the Kaplan-Meier method
## applied to the censored times reversing the roles of event status and censored.
##
## Table of quantiles and corresponding confidence limits:
## q quantile lower upper
## <num> <num> <num> <num>
## 1: 0.00 NA NA NA
## 2: 0.25 9.1 8.4 9.8
## 3: 0.50 6.5 6.1 7.0
## 4: 0.75 4.4 3.9 4.8
## 5: 1.00 1.5 1.5 2.0
##
##
## Median with interquartile range (IQR):
## Median (IQR)
## <char>
## 1: 6.48 (4.36;9.14)
The restricted mean survival time (RMST) is the area under the Kaplan-Meier curve up to a chosen time point. It is the average number of years a patient is alive during that period, for example during the first 10 years. For a quick guide, see Han & Jung (2022).
A hazard ratio is a relative measure that assumes the effect is constant over time (proportional hazards) and is hard to translate into an absolute effect. The RMST is in the unit of time (years gained or lost), so it is easy to interpret, it does not need the proportional hazards assumption, and it can still show a difference when the survival curves cross. The difference in RMST between two groups (the area between their Kaplan-Meier curves) also gives a better picture of the absolute treatment benefit than the risk difference at a single time point, which can lead to treating the wrong patients. For that reason it is recommended for assessing absolute treatment effects in time-to-event data (Maas et al., 2025).
You can get it from the summary() of the Kaplan-Meier
fit by setting rmean to the time point. The result contains
a table with the RMST (rmean) and its standard
error (se(rmean)) per group. We use the outcome restricted
to the 10-year horizon:
# restricted mean survival time up to the horizon (10 years)
# $table takes the table with the RMST out of the summary
summary(KM_10, rmean=horizon)$table
## records n.max n.start events rmean se(rmean) median
## trt=D-penicillamine 158 158 158 63 7.148479 0.2827057 8.991781
## trt=Placebo 154 154 154 57 7.285271 0.2954189 9.391781
## 0.95LCL 0.95UCL
## trt=D-penicillamine 7.076712 NA
## trt=Placebo 8.465753 NA
# extract the RMST (and its standard error) for one treatment arm:
# [row name, column names] picks one row and two columns of the table
summary(KM_10, rmean=horizon)$table["trt=Placebo", c("rmean", "se(rmean)")]
## rmean se(rmean)
## 7.2852709 0.2954189
For example, an rmean of 7.1 means that patients in that
group live on average 7.1 of the first 10 years. The RMST needs a time
point to restrict at, otherwise R restricts at the largest observed
time. This is another reason to limit the outcome to the horizon.
Try each exercise yourself first. Then open the field below it to see the answer, the code and its results.
Exercise 1. What are the Kaplan-Meier estimates at 5 years? How do they compare with the estimates at the 10-year horizon?
A Kaplan-Meier estimate at 5 years does not depend on follow-up after
5 years, so you do not need to define a new outcome. Use
summary() of the fit with S_10 and set
time=5:
# Kaplan-Meier estimates at 5 years
summary(KM_10, time=5)
## Call: survfit(formula = S_10 ~ trt, data = data)
##
## 106 observations deleted due to missingness
## trt=D-penicillamine
## time n.risk n.event survival std.err lower 95% CI
## 5.0000 82.0000 43.0000 0.7077 0.0379 0.6371
## upper 95% CI
## 0.7861
##
## trt=Placebo
## time n.risk n.event survival std.err lower 95% CI
## 5.0000 77.0000 42.0000 0.7146 0.0376 0.6445
## upper 95% CI
## 0.7923
Compare these estimates with the 10-year estimates above. A Kaplan-Meier curve can only go down over time, so the survival estimate at 5 years is at least as high as the estimate at 10 years.
Exercise 2. Calculate the naive median,
median(data$time)/365, and compare it with the reverse
Kaplan-Meier estimate. Which one would you report, and why?
# naive median of the observed times (in years)
median(data$time)/365
## [1] 4.739726
# median follow-up with the reverse Kaplan-Meier method (in years)
stats::quantile(reverse_KM, probs=0.5)
## $quantile
## 50
## 6.479452
##
## $lower
## 50
## 6.093151
##
## $upper
## 50
## 7.021918
Report the reverse Kaplan-Meier estimate. The naive median mixes patients who were followed until the end of the study with patients whose follow-up stopped early because they died, which makes the typical follow-up look shorter than it is.
Exercise 3. Calculate the restricted mean survival time (RMST) up to 5 years instead of 10 years for each treatment group, and extract the RMST of the placebo arm only. How many of the first 5 years do patients on average live?
# RMST up to 5 years for both treatment arms
summary(KM_10, rmean=5)$table
## records n.max n.start events rmean se(rmean) median
## trt=D-penicillamine 158 158 158 63 4.30216 0.1060121 8.991781
## trt=Placebo 154 154 154 57 4.18246 0.1190897 9.391781
## 0.95LCL 0.95UCL
## trt=D-penicillamine 7.076712 NA
## trt=Placebo 8.465753 NA
# extract the RMST (and its standard error) of the placebo arm only
summary(KM_10, rmean=5)$table["trt=Placebo", c("rmean", "se(rmean)")]
## rmean se(rmean)
## 4.1824596 0.1190897
The fit KM_10 covers 10 years, but rmean=5
restricts the area under the curve to the first 5 years, so you do not
need a new fit. The RMST up to 5 years is at most 5, and smaller than
the RMST up to 10 years.
Exercise 4. Why do we truncate the outcome at the horizon instead of just using all follow-up?
# how many patients are followed for longer than the 10-year horizon?
table(S[, 1] > horizon)
##
## FALSE TRUE
## 383 35
A prediction model gives the risk at a specific time point, here 10 years, and the baseline hazard is estimated at that time point. Later follow-up is based on few patients at risk and would not be part of the risk that we predict. Truncating at the horizon keeps the outcome consistent with the prediction.
Most clinical data sets have missing values. Deleting every patient with a missing value (complete case analysis) wastes data and can bias the model. Multiple imputation fills in missing values several times, using the observed data, so that the uncertainty due to missingness is reflected in the results. In this part you visualize the missing data, build an imputation model that suits survival data, and check the imputed values.
Before imputing, look at how much is missing and where.
# visualize missingness: every column is a variable, every row a patient,
# and the black cells are missing values
visdat::vis_miss(data)
The imputation model must contain the outcome, otherwise the imputed
predictors are unrelated to the outcome and the model’s effects are
biased towards zero. For time-to-event data, it is not
a good idea to put time and status in the
imputation model directly. The recommended approach (White &
Royston, 2009) is to include:
We calculate both on the same outcome that we model later: death
within the 10-year horizon (S_10), with transplant treated
as censored. The raw time and status variables
are then left out of the imputation model.
# outcome information for the imputation model, based on the analysis outcome S_10
data_mi <- data # work on a copy, so that `data` stays as it is
data_mi$time <- NULL # NULL removes a column
data_mi$status <- NULL
# 1. the event indicator (0/1) of the outcome within 10 years
data_mi$event_10 <- S_10[, "status"]
# 2. the Nelson-Aalen estimate of the cumulative hazard at each patient's time
data_mi$cum_hazard <- mice::nelsonaalen(data.frame(time=S_10[, "time"],
status=S_10[, "status"]),
time, status)
We create m = 5 imputed data sets. seed=1
makes the result reproducible.
Note: five imputations is on the small side. A rule of thumb is to use at least as many imputations as the percentage of incomplete cases (White et al., 2011). Let’s check that percentage here:
# percentage of patients with at least one missing value
round(100*mean(!complete.cases(data_mi)), 1)
## [1] 34
With this percentage, m should probably be larger than 5
(maybe 40 or more). For this tutorial we keep m = 5 so that
all following steps run quickly, but in your own project you should use
a larger m based on the fraction of incomplete cases.
# multiple imputation (m should maybe be larger, see the text above)
m <- 5 # number of imputed data sets
# mice() imputes the missing values: m=number of data sets, seed=1 makes the
# result reproducible, print=FALSE hides the progress messages
imputed_data <- mice::mice(data_mi, m=m, seed=1, print=FALSE)
# problems that mice ran into (NULL means none)
imputed_data$loggedEvents
## NULL
loggedEvents lists problems that mice ran
into; NULL means none.
The fraction of missing information (FMI) is the part of the
uncertainty (variance) of an estimate that is due to missing data. It
can be much lower or higher than the percentage of missing values,
because it depends on how well the other variables predict the missing
ones. We calculate it per variable: fit an intercept-only model on every
imputed data set (its estimate is the mean, or the log odds for a binary
variable), pool the results with Rubin’s rules using
pool(), and take the FMI from the pooled result. We skip
categorical variables with more than 2 levels, because a single
intercept does not describe them.
To do this for every variable we need two basic tools: a function and a for loop.
A function is a piece of code that you give a name, so that you can run it again and again with different input. For example, a function that squares a number:
# define a function: x is the input, the last line is the output
square <- function(x){
x^2
}
square(3)
## [1] 9
square(10)
## [1] 100
A for loop repeats code for every element of a vector. For example, print the square of each number in a vector:
# the loop runs once for every value of `number`
for (number in c(2, 3, 4)){
print(square(number))
}
## [1] 4
## [1] 9
## [1] 16
Now we write a function that calculates the FMI for one variable, and use a for loop to run it for every variable with missing values.
# ---- 1. Define the function -------------------------------------------------
# A function is defined with function(argument){ ... }. Here the argument `var`
# is the name of one variable (for example "chol"), and the function returns a
# small table with the percentage missing and the FMI of that variable.
calculate_fmi <- function(var){
# data_mi[[var]] is the column with the name stored in `var`.
# is.na() gives TRUE for missing values, mean() of TRUE/FALSE gives the
# proportion of TRUE, and 100* turns it into a percentage (rounded to 1 decimal)
pct_missing <- round(100*mean(is.na(data_mi[[var]])), 1)
# We fit a model with only an intercept ("var ~ 1") on each of the m imputed
# data sets. with(imputed_data, ...) runs the model on every imputed data set.
# paste(var, "~ 1") builds the formula as text, for example "chol ~ 1", and
# as.formula() turns the text into a formula.
if (is.factor(data_mi[[var]]) && nlevels(data_mi[[var]]) == 2){
# binary variable: logistic regression (the intercept is the log odds)
fit <- with(imputed_data, glm(as.formula(paste(var, "~ 1")), family=binomial))
} else if (is.numeric(data_mi[[var]])){
# continuous variable: linear regression (the intercept is the mean)
fit <- with(imputed_data, lm(as.formula(paste(var, "~ 1"))))
} else {
# other variables (categorical with more than 2 levels): no FMI, so we
# return a row with NA and stop here (return() ends the function)
return(data.frame(variable=var, pct_missing=pct_missing, fmi=NA))
}
# pool() combines the m models into one estimate with Rubin's rules. The
# pooled result contains the fraction of missing information (`fmi`), a
# proportion between 0 and 1; we show it as a percentage.
pooled <- mice::pool(fit)
# The last line of a function is its result: a one-row table for this variable
data.frame(variable=var, pct_missing=pct_missing,
fmi=round(100*pooled$pooled$fmi, 1))
}
# ---- 2. Use the function in a for loop ----------------------------------------
# colSums(is.na(data_mi)) counts the missing values per column; we keep the names
# of the columns that have at least one missing value.
vars_missing <- names(data_mi)[colSums(is.na(data_mi)) > 0]
# Start with an empty table and add one row per variable. In every round of the
# loop `var` is the next name in `vars_missing`; rbind() puts the new row under
# the rows that we already have.
fmi_results <- data.frame()
for (var in vars_missing){
fmi_results <- rbind(fmi_results, calculate_fmi(var))
}
# print the table
fmi_results
A high FMI means that the missing data add a lot of uncertainty and
that more imputations are needed to estimate it stably. A common rule of
thumb is that m should be at least 100 times the largest
FMI (as a proportion, so an FMI of 0.3 gives m = 30) (White
et al., 2011).
Compare the distribution of observed values (blue) with the imputed values (red) for each imputation. They do not need to be identical (missingness may depend on other variables), but strange differences deserve a closer look.
# visualize imputations
mice::densityplot(imputed_data)
Try each exercise yourself first. Then open the field below it to see the answer, the code and its results.
Exercise 1. Which variables have the most missing values, and what percentage of patients have at least one missing value?
# percentage missing per variable, most missing first
round(100*sort(colMeans(is.na(data)), decreasing=TRUE), 1)
## trig chol copper trt ascites hepato spiders alk.phos
## 32.5 32.1 25.8 25.4 25.4 25.4 25.4 25.4
## ast platelet stage protime time status age sex
## 25.4 2.6 1.4 0.5 0.0 0.0 0.0 0.0
## edema bili albumin center
## 0.0 0.0 0.0 0.0
# percentage of patients with at least one missing value
round(100*mean(!complete.cases(data_mi)), 1)
## [1] 34
The first table gives the percentage of missing values per variable, the second number the percentage of patients with at least one missing value. The second number is the one used for the rule of thumb for the number of imputations.
Exercise 2. Look at the density plots. Do the
imputed values (red) look plausible compared with the observed values
(blue) for chol, trig and
copper?
mice::densityplot(imputed_data, ~ chol)
mice::densityplot(imputed_data, ~ trig)
mice::densityplot(imputed_data, ~ copper)
Compare the shape of the red lines (imputed) with the blue line
(observed). These variables are right-skewed, and the imputed values do
not need to match the observed values exactly: if missingness depends on
other variables, such as being in the trial or not, the imputed values
may legitimately differ. By default mice uses predictive
mean matching for continuous variables, which only imputes values that
occur in the observed data, so impossible values (for example negative
cholesterol) cannot appear. Be suspicious if the imputed distribution is
completely different or has a strange shape.
Exercise 3. Increase m to the
rule-of-thumb value and re-run. How long does it take, and do the
density plots change?
This can take a minute or two. The new imputations are stored in
imputed_data_big, so the tutorial keeps using
m = 5.
# rule of thumb: at least as many imputations as the percentage of incomplete cases
m_rule <- ceiling(100*mean(!complete.cases(data_mi)))
m_rule
## [1] 34
# impute again with more imputations and time it
time_big <- system.time(
imputed_data_big <- mice::mice(data_mi, m=m_rule, seed=1, print=FALSE)
)
time_big["elapsed"]
## elapsed
## 18.85
# density plot with more imputations
mice::densityplot(imputed_data_big, ~ chol)
The time increases roughly in proportion to the number of imputations. The density plots look similar (more red lines), because the imputation model is the same. The benefit of more imputations is more stable estimates and standard errors of the final model, not a different imputation distribution.
Exercise 4. Suppose the outcome is missing for some
patients (for example, we lost track of their vital status). Should you
delete these patients before imputation? Show how to impute with the
outcome in the imputation model and delete the patients with a missing
outcome afterwards. To mimic this, set event_10 to missing
for every 20th patient.
Do not delete them before imputation: these patients still have observed predictor values that help to impute the missing predictors of other patients. Instead (1) first create a variable that flags who has a missing outcome, (2) impute with the outcome in the imputation model, and (3) delete the flagged rows from every imputed data set afterwards. This is called multiple imputation, then deletion (von Hippel, 2007).
# mimic missing outcomes: set event_10 to missing for every 20th patient
data_mi_out <- data_mi
data_mi_out$event_10[seq(1, nrow(data_mi_out), by=20)] <- NA
# 1. flag the patients with a missing outcome BEFORE imputation
data_mi_out$outcome_missing <- is.na(data_mi_out$event_10)
# 2. impute, with the outcome in the imputation model
# (the flag itself is not used to impute other variables)
pred_matrix <- mice::make.predictorMatrix(data_mi_out)
pred_matrix[, "outcome_missing"] <- 0
imputed_out <- mice::mice(data_mi_out, m=5, predictorMatrix=pred_matrix,
seed=1, print=FALSE)
# 3. after imputation: delete the rows with a missing outcome in each imputed data set
imputed_sets <- mice::complete(imputed_out, "all")
for (i in seq_along(imputed_sets)){
imputed_sets[[i]] <- imputed_sets[[i]][!imputed_sets[[i]]$outcome_missing, ]
}
sapply(imputed_sets, nrow)
## 1 2 3 4 5
## 397 397 397 397 397
Each imputed data set now has the patients with an observed outcome only, and the predictors of these patients were imputed using all patients. The imputed outcomes themselves are not analysed. This approach is for a missing outcome only; always keep the outcome in the imputation model.
In this part you fit a Cox proportional hazards model on the imputed
data, check whether continuous predictors need a non-linear effect
(using rcs()), select predictors with backward selection,
and present the final model (effects, variable importance and a
nomogram). cph (from rms) fits a Cox model
with extra tools for prediction models, and fit.mult.impute
fits the model on every imputed data set and pools the results with
Rubin’s rules.
datadist stores the ranges of the predictors, which
rms needs for summaries, plots and nomograms.
# set data distributions: datadist() stores the ranges of all predictors, and
# options() tells rms to use it (it looks for the object called "dd")
dd <- rms::datadist(data)
options(datadist="dd")
rcs()Continuous predictors such as age, bilirubin (bili) and
albumin are often assumed to have a linear effect: every extra
year or unit changes the log hazard by the same amount. That is rarely
exactly true. A flexible way to model a non-linear effect is a
restricted cubic spline, rcs(x, k) from
rms:
k is the number of knots, the points where the
curve is allowed to bend. With 3 knots, they are placed at the 10th,
50th and 90th percentile of the variable (with more knots, at other
percentiles). Use 3 to 5 knots: more knots make the curve more flexible,
but need more parameters and can overfit.rcs(x, 3) uses 2 coefficients: one for x
itself (the linear part) and one for x' (the non-linear
part). rcs(x, 4) uses 3, and so on. A linear effect uses
1.We first model all continuous predictors flexibly with
rcs(..., 3) and test whether the non-linear part is
significant. Fitting flexible effects for all continuous variables is a
good way to see if linearity is a reasonable assumption.
# fit full model with all candidate variables on each imputed data set
# model all continuous variables flexible to see if that is needed
# fit.mult.impute() arguments:
# - the formula: outcome ~ predictor 1 + predictor 2 + ...
# - cph: the function that fits the Cox model on each imputed data set
# - fitargs: extra arguments for cph (keep x, y and the survival estimates)
# - xtrans: the imputed data sets
# - data: the original data, with the missing values
# - fit.reps=TRUE: keep the model of every imputed data set
full_non_linear_model <- Hmisc::fit.mult.impute(S_10 ~ rcs(age, 3) +
rcs(bili, 3) +
rcs(albumin, 3) + sex +
ascites + hepato + stage,
cph,
fitargs=list(x=TRUE, y=TRUE, surv=TRUE),
xtrans=imputed_data,
data=data,
fit.reps=TRUE)
##
## Wald Statistic Information
##
## Variance Inflation Factors Due to Imputation:
##
## age age' bili bili' albumin albumin'
## 1.01 1.01 1.01 1.00 1.00 1.01
## sex=Female ascites=Yes hepato=Yes stage=2 stage=3 stage=4
## 1.00 1.15 1.09 1.09 1.08 1.09
##
## Fraction of Missing Information:
##
## age age' bili bili' albumin albumin'
## 0.01 0.01 0.01 0.00 0.00 0.01
## sex=Female ascites=Yes hepato=Yes stage=2 stage=3 stage=4
## 0.00 0.13 0.08 0.08 0.07 0.08
##
## d.f. for t-distribution for Tests of Single Coefficients:
##
## age age' bili bili' albumin albumin'
## 70992.09 29693.89 106704.51 199948.55 205552.84 32084.48
## sex=Female ascites=Yes hepato=Yes stage=2 stage=3 stage=4
## 358897.26 242.90 630.61 651.50 799.20 602.05
##
## The following fit components were averaged over the 5 model fits:
##
## linear.predictors means stats center
summary(full_non_linear_model)
## Effects Response : S_10
##
## Factor Low High Diff. Effect S.E. Lower 0.95
## age 42.8320 58.241 15.4090 0.556350 0.15466 0.2532300
## Hazard Ratio 42.8320 58.241 15.4090 1.744300 NA 1.2882000
## bili 0.8000 3.400 2.6000 1.495800 0.20069 1.1025000
## Hazard Ratio 0.8000 3.400 2.6000 4.462900 NA 3.0116000
## albumin 3.2425 3.770 0.5275 -0.360650 0.12539 -0.6064000
## Hazard Ratio 3.2425 3.770 0.5275 0.697220 NA 0.5453100
## sex - Male:Female 2.0000 1.000 NA 0.081707 0.25309 -0.4143400
## Hazard Ratio 2.0000 1.000 NA 1.085100 NA 0.6607700
## ascites - Yes:No 1.0000 2.000 NA 0.438920 0.26828 -0.0868970
## Hazard Ratio 1.0000 2.000 NA 1.551000 NA 0.9167700
## hepato - No:Yes 2.0000 1.000 NA -0.087502 0.22035 -0.5193800
## Hazard Ratio 2.0000 1.000 NA 0.916220 NA 0.5948900
## stage - 1:3 3.0000 1.000 NA -0.730420 0.73833 -2.1775000
## Hazard Ratio 3.0000 1.000 NA 0.481710 NA 0.1133200
## stage - 2:3 3.0000 2.000 NA -0.482070 0.27063 -1.0125000
## Hazard Ratio 3.0000 2.000 NA 0.617500 NA 0.3633100
## stage - 4:3 3.0000 4.000 NA 0.415320 0.21663 -0.0092709
## Hazard Ratio 3.0000 4.000 NA 1.514900 NA 0.9907700
## Upper 0.95
## 0.859470
## 2.361900
## 1.889100
## 6.613600
## -0.114890
## 0.891460
## 0.577760
## 1.782000
## 0.964730
## 2.624100
## 0.344380
## 1.411100
## 0.716690
## 2.047600
## 0.048357
## 1.049500
## 0.839920
## 2.316200
The anova table displays chi-square tests for
non-linearity. Look at the row “Nonlinear” under each continuous
predictor: a small p-value means that the effect is not a straight
line.
# display chi-square tests for non-linearities
anova(full_non_linear_model)
## Wald Statistics Response: S_10
##
## Factor Chi-Square d.f. P
## age 13.11 2 0.0014
## Nonlinear 3.51 1 0.0609
## bili 95.92 2 <.0001
## Nonlinear 35.10 1 <.0001
## albumin 11.80 2 0.0027
## Nonlinear 0.12 1 0.7263
## sex 0.10 1 0.7468
## ascites 2.68 1 0.1018
## hepato 0.16 1 0.6913
## stage 11.06 3 0.0114
## TOTAL NONLINEAR 36.64 3 <.0001
## TOTAL 215.30 12 <.0001
Here bilirubin has a clearly non-linear effect, age is borderline, and albumin is close to linear. Plot the estimated effects to see what this looks like:
# display predictor effects (log hazard) for the flexible model
plot(Predict(full_non_linear_model))
The curve for bilirubin rises steeply at low values and flattens out
at high values, which a straight line would miss. Based on these results
we keep rcs() for age and bilirubin and a linear effect for
albumin in the full model.
Now we fit the full model with rcs() for age and
bilirubin, a linear effect for albumin, and the other candidate
predictors. This model is the starting point for variable selection.
# fit full model with all candidate variables on each imputed data set
# we store the formula in an object, so that we can reuse it later (Tutorial 2)
# rcs(age, 3): spline for age; albumin without rcs(): a straight line
form_full <- S_10 ~ rcs(age, 3) + rcs(bili, 3) + albumin + sex +
ascites + hepato + stage
full_model <- Hmisc::fit.mult.impute(form_full,
cph,
fitargs=list(x=TRUE, y=TRUE, surv=TRUE),
xtrans=imputed_data,
data=data,
fit.reps=TRUE)
##
## Wald Statistic Information
##
## Variance Inflation Factors Due to Imputation:
##
## age age' bili bili' albumin sex=Female
## 1.01 1.01 1.01 1.00 1.00 1.00
## ascites=Yes hepato=Yes stage=2 stage=3 stage=4
## 1.14 1.10 1.08 1.08 1.09
##
## Fraction of Missing Information:
##
## age age' bili bili' albumin sex=Female
## 0.01 0.01 0.01 0.00 0.00 0.00
## ascites=Yes hepato=Yes stage=2 stage=3 stage=4
## 0.13 0.09 0.08 0.07 0.08
##
## d.f. for t-distribution for Tests of Single Coefficients:
##
## age age' bili bili' albumin sex=Female
## 66673.51 25006.17 103533.09 203798.25 31013743.69 337198.38
## ascites=Yes hepato=Yes stage=2 stage=3 stage=4
## 251.90 492.12 655.56 784.71 576.02
##
## The following fit components were averaged over the 5 model fits:
##
## linear.predictors means stats center
# display predictor effects (hazard ratios, comparing the 25th and 75th percentile)
summary(full_model)
## Effects Response : S_10
##
## Factor Low High Diff. Effect S.E. Lower 0.95
## age 42.8320 58.241 15.4090 0.557540 0.15483 0.254070
## Hazard Ratio 42.8320 58.241 15.4090 1.746400 NA 1.289300
## bili 0.8000 3.400 2.6000 1.502000 0.20048 1.109100
## Hazard Ratio 0.8000 3.400 2.6000 4.490800 NA 3.031600
## albumin 3.2425 3.770 0.5275 -0.380950 0.11312 -0.602660
## Hazard Ratio 3.2425 3.770 0.5275 0.683210 NA 0.547360
## sex - Male:Female 2.0000 1.000 NA 0.084776 0.25332 -0.411730
## Hazard Ratio 2.0000 1.000 NA 1.088500 NA 0.662500
## ascites - Yes:No 1.0000 2.000 NA 0.453960 0.26293 -0.061373
## Hazard Ratio 1.0000 2.000 NA 1.574500 NA 0.940470
## hepato - No:Yes 2.0000 1.000 NA -0.081244 0.22032 -0.513070
## Hazard Ratio 2.0000 1.000 NA 0.921970 NA 0.598660
## stage - 1:3 3.0000 1.000 NA -0.731010 0.73871 -2.178800
## Hazard Ratio 3.0000 1.000 NA 0.481420 NA 0.113170
## stage - 2:3 3.0000 2.000 NA -0.483820 0.27056 -1.014100
## Hazard Ratio 3.0000 2.000 NA 0.616420 NA 0.362720
## stage - 4:3 3.0000 4.000 NA 0.412780 0.21622 -0.010994
## Hazard Ratio 3.0000 4.000 NA 1.511000 NA 0.989070
## Upper 0.95
## 0.86101
## 2.36550
## 1.89500
## 6.65230
## -0.15924
## 0.85279
## 0.58128
## 1.78830
## 0.96929
## 2.63610
## 0.35058
## 1.41990
## 0.71683
## 2.04790
## 0.04647
## 1.04760
## 0.83656
## 2.30840
# display chi-square tests for variable importance: a larger chi-square
# (relative to its degrees of freedom) means a more important predictor
anova(full_model)
## Wald Statistics Response: S_10
##
## Factor Chi-Square d.f. P
## age 13.18 2 0.0014
## Nonlinear 3.70 1 0.0546
## bili 95.81 2 <.0001
## Nonlinear 35.72 1 <.0001
## albumin 11.34 1 0.0008
## sex 0.11 1 0.7379
## ascites 2.98 1 0.0842
## hepato 0.14 1 0.7123
## stage 11.03 3 0.0116
## TOTAL NONLINEAR 36.34 2 <.0001
## TOTAL 214.11 11 <.0001
We check the proportional hazards assumption, which says that the effect of a predictor (the hazard ratio) is constant over time. A small p-value for a predictor suggests a violation.
# test proportional hazards assumption: a small p-value suggests that the
# effect of that predictor changes over time
survival::cox.zph(full_model)
## chisq df p
## rcs(age, 3) 0.5573 66674 1
## rcs(bili, 3) 0.4378 25006 1
## albumin 0.6281 103533 1
## sex 0.5998 203798 1
## ascites 2.2708 31013744 1
## hepato 0.0879 337198 1
## stage 4.3057 252 1
## GLOBAL 6.8567 31752713 1
# display predictor effects: Predict() calculates the effect of each predictor
# on the log hazard (others held at their median), plot() draws it
plot(Predict(full_model))
We use backward selection with fastbw(): starting from
the full model, variables with a p-value above 0.05 are removed one by
one. type="individual" tests each variable separately. Note
that rcs(age, 3) is one variable with two terms
(age and age'), so it is listed once, as
age.
For this illustrative exercise fastbw() is fine. There
are better ways to select variables, for example penalized regression
(LASSO or elastic net) or relying on clinical knowledge instead of
p-values, as backward selection tends to overfit. This is also why we
estimate a shrinkage factor in Tutorial 2.
# perform backward selection with fastbw():
# type="individual": test each variable separately
# rule="p": remove variables based on their p-value
# sls=0.05: remove a variable if its p-value is above 0.05
fastbw(full_model, type="individual", rule="p", sls=0.05)
##
## Deleted Chi-Sq d.f. P Residual d.f. P AIC
## sex 0.11 1 0.7379 0.11 1 0.7379 -1.89
## hepato 0.15 1 0.6988 0.26 2 0.8773 -3.74
## ascites 3.05 1 0.0807 3.31 3 0.3460 -2.69
##
## Approximate Estimates after Deleting Factors
##
## Coef S.E. Wald Z P
## age 0.06770 0.02400 2.8207 4.792e-03
## age' -0.03661 0.02257 -1.6224 1.047e-01
## bili 0.79304 0.10919 7.2627 3.796e-13
## bili' -2.30232 0.37783 -6.0936 1.104e-09
## albumin -0.83765 0.20184 -4.1501 3.323e-05
## stage=2 0.29603 0.74623 0.3967 6.916e-01
## stage=3 0.73805 0.72262 1.0214 3.071e-01
## stage=4 1.23083 0.72734 1.6922 9.060e-02
##
## Factors in Final Model
##
## [1] age bili albumin stage
# this gives final model: age, bili, albumin, stage
We then fit the backward selected model.
# fit backward selected model: the same as before, with only the selected variables
form_bw <- S_10 ~ rcs(age, 3) + rcs(bili, 3) + albumin + stage
bw_model <- Hmisc::fit.mult.impute(form_bw,
cph,
fitargs=list(x=TRUE, y=TRUE, surv=TRUE),
xtrans=imputed_data,
data=data,
fit.reps=TRUE)
##
## Wald Statistic Information
##
## Variance Inflation Factors Due to Imputation:
##
## age age' bili bili' albumin stage=2 stage=3 stage=4
## 1.00 1.00 1.00 1.00 1.01 1.09 1.06 1.07
##
## Fraction of Missing Information:
##
## age age' bili bili' albumin stage=2 stage=3 stage=4
## 0.00 0.00 0.00 0.00 0.01 0.08 0.06 0.06
##
## d.f. for t-distribution for Tests of Single Coefficients:
##
## age age' bili bili' albumin stage=2
## 20946309.50 131064158.78 381202.54 364582.19 97205.48 592.59
## stage=3 stage=4
## 1199.64 1069.91
##
## The following fit components were averaged over the 5 model fits:
##
## linear.predictors means stats center
# display predictor effects
summary(bw_model)
## Effects Response : S_10
##
## Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
## age 42.8320 58.241 15.4090 0.59769 0.14878 0.306080 0.88930
## Hazard Ratio 42.8320 58.241 15.4090 1.81790 NA 1.358100 2.43340
## bili 0.8000 3.400 2.6000 1.51860 0.19498 1.136400 1.90070
## Hazard Ratio 0.8000 3.400 2.6000 4.56560 NA 3.115500 6.69060
## albumin 3.2425 3.770 0.5275 -0.44221 0.10613 -0.650230 -0.23419
## Hazard Ratio 3.2425 3.770 0.5275 0.64261 NA 0.521920 0.79121
## stage - 1:3 3.0000 1.000 NA -0.71732 0.72234 -2.133100 0.69845
## Hazard Ratio 3.0000 1.000 NA 0.48806 NA 0.118470 2.01060
## stage - 2:3 3.0000 2.000 NA -0.45077 0.27036 -0.980680 0.07913
## Hazard Ratio 3.0000 2.000 NA 0.63713 NA 0.375060 1.08230
## stage - 4:3 3.0000 4.000 NA 0.48047 0.19963 0.089196 0.87175
## Hazard Ratio 3.0000 4.000 NA 1.61680 NA 1.093300 2.39110
The effect of age per decade. Because age has a non-linear effect (a spline), the hazard ratio for 10 years of age depends on where you start. We therefore compare ages within the observed range, here 50 versus 60 years. (A comparison of age 0 with age 10 would be outside the range of the data and meaningless.)
# present effect for age per decade (from 50 to 60 years)
# age=c(50, 60) sets the two values that are compared
summary(bw_model, age=c(50, 60))
## Effects Response : S_10
##
## Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
## age 50.0000 60.00 10.0000 0.24457 0.094077 0.060179 0.42895
## Hazard Ratio 50.0000 60.00 10.0000 1.27710 NA 1.062000 1.53560
## bili 0.8000 3.40 2.6000 1.51860 0.194980 1.136400 1.90070
## Hazard Ratio 0.8000 3.40 2.6000 4.56560 NA 3.115500 6.69060
## albumin 3.2425 3.77 0.5275 -0.44221 0.106130 -0.650230 -0.23419
## Hazard Ratio 3.2425 3.77 0.5275 0.64261 NA 0.521920 0.79121
## stage - 1:3 3.0000 1.00 NA -0.71732 0.722340 -2.133100 0.69845
## Hazard Ratio 3.0000 1.00 NA 0.48806 NA 0.118470 2.01060
## stage - 2:3 3.0000 2.00 NA -0.45077 0.270360 -0.980680 0.07913
## Hazard Ratio 3.0000 2.00 NA 0.63713 NA 0.375060 1.08230
## stage - 4:3 3.0000 4.00 NA 0.48047 0.199630 0.089196 0.87175
## Hazard Ratio 3.0000 4.00 NA 1.61680 NA 1.093300 2.39110
Variable importance with the chi-square test statistic (minus its degrees of freedom is a fair way to compare):
# display variable importance using chi-square test statistic
anova(bw_model)
## Wald Statistics Response: S_10
##
## Factor Chi-Square d.f. P
## age 16.27 2 0.0003
## Nonlinear 2.65 1 0.1033
## bili 116.25 2 <.0001
## Nonlinear 37.53 1 <.0001
## albumin 17.36 1 <.0001
## stage 14.66 3 0.0021
## TOTAL NONLINEAR 37.70 2 <.0001
## TOTAL 206.90 8 <.0001
# visualize predictor effects
plot(Predict(bw_model))
A nomogram lets you calculate the risk of an individual patient by hand.
# nomogram: nomogram() calculates the points of every predictor value,
# plot() draws it
plot(rms::nomogram(bw_model))
Try each exercise yourself first. Then open the field below it to see the answer, the code and its results. The answers use new object names, so the saved models do not change.
Exercise 1. In the nomogram, what is the total number of points for a 75-year-old patient with a bilirubin of 2, an albumin of 3.6 and stage 4? Can you tell which variable contributes the most?
Printing the nomogram gives the points for every value of every predictor (for age, in steps of 5 years):
nom <- rms::nomogram(bw_model)
print(nom)
## Points per unit of linear predictor: 25.99955
## Linear predictor units per point : 0.03846221
##
##
## age Points
## 25 0
## 30 9
## 35 18
## 40 26
## 45 34
## 50 41
## 55 45
## 60 47
## 65 49
## 70 50
## 75 51
## 80 52
##
##
## bili Points
## 0 0
## 2 38
## 4 61
## 6 71
## 8 75
## 10 77
## 12 80
## 14 82
## 16 85
## 18 87
## 20 90
## 22 92
## 24 95
## 26 97
## 28 100
##
##
## albumin Points
## 1.8 65
## 2.0 61
## 2.2 57
## 2.4 52
## 2.6 48
## 2.8 44
## 3.0 39
## 3.2 35
## 3.4 31
## 3.6 26
## 3.8 22
## 4.0 17
## 4.2 13
## 4.4 9
## 4.6 4
## 4.8 0
##
##
## stage Points
## 1 0
## 2 7
## 3 19
## 4 31
Look up the points for age 75, bilirubin 2, albumin 3.6 and stage 4, and add them up. The total is read on the “Total Points” axis, which converts the points to a risk. The variable with the widest range of points (the longest axis in the nomogram) contributes the most to the prediction.
Exercise 2. Change the age comparison to
c(40, 50) and c(60, 70). Why do the hazard
ratios differ?
summary(bw_model, age=c(40, 50))
## Effects Response : S_10
##
## Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
## age 40.0000 50.00 10.0000 0.55853 0.17340 0.218660 0.89839
## Hazard Ratio 40.0000 50.00 10.0000 1.74810 NA 1.244400 2.45570
## bili 0.8000 3.40 2.6000 1.51860 0.19498 1.136400 1.90070
## Hazard Ratio 0.8000 3.40 2.6000 4.56560 NA 3.115500 6.69060
## albumin 3.2425 3.77 0.5275 -0.44221 0.10613 -0.650230 -0.23419
## Hazard Ratio 3.2425 3.77 0.5275 0.64261 NA 0.521920 0.79121
## stage - 1:3 3.0000 1.00 NA -0.71732 0.72234 -2.133100 0.69845
## Hazard Ratio 3.0000 1.00 NA 0.48806 NA 0.118470 2.01060
## stage - 2:3 3.0000 2.00 NA -0.45077 0.27036 -0.980680 0.07913
## Hazard Ratio 3.0000 2.00 NA 0.63713 NA 0.375060 1.08230
## stage - 4:3 3.0000 4.00 NA 0.48047 0.19963 0.089196 0.87175
## Hazard Ratio 3.0000 4.00 NA 1.61680 NA 1.093300 2.39110
summary(bw_model, age=c(60, 70))
## Effects Response : S_10
##
## Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
## age 60.0000 70.00 10.0000 0.10094 0.15461 -0.202090 0.40397
## Hazard Ratio 60.0000 70.00 10.0000 1.10620 NA 0.817020 1.49780
## bili 0.8000 3.40 2.6000 1.51860 0.19498 1.136400 1.90070
## Hazard Ratio 0.8000 3.40 2.6000 4.56560 NA 3.115500 6.69060
## albumin 3.2425 3.77 0.5275 -0.44221 0.10613 -0.650230 -0.23419
## Hazard Ratio 3.2425 3.77 0.5275 0.64261 NA 0.521920 0.79121
## stage - 1:3 3.0000 1.00 NA -0.71732 0.72234 -2.133100 0.69845
## Hazard Ratio 3.0000 1.00 NA 0.48806 NA 0.118470 2.01060
## stage - 2:3 3.0000 2.00 NA -0.45077 0.27036 -0.980680 0.07913
## Hazard Ratio 3.0000 2.00 NA 0.63713 NA 0.375060 1.08230
## stage - 4:3 3.0000 4.00 NA 0.48047 0.19963 0.089196 0.87175
## Hazard Ratio 3.0000 4.00 NA 1.61680 NA 1.093300 2.39110
# plot the effect of age on the log hazard
plot(Predict(bw_model, age))
The plot shows why: age has a non-linear effect (modelled with
rcs()), so the effect of age is not a straight line. The
change in log hazard for a 10-year difference therefore depends on where
you start. With a linear effect, the hazard ratio per decade would be
the same everywhere.
Exercise 3. Change sls=0.05 to
sls=0.157 (similar to AIC). Does the selected model
change?
fastbw(full_model, type="individual", rule="p", sls=0.157)
##
## Deleted Chi-Sq d.f. P Residual d.f. P AIC
## sex 0.11 1 0.7379 0.11 1 0.7379 -1.89
## hepato 0.15 1 0.6988 0.26 2 0.8773 -3.74
##
## Approximate Estimates after Deleting Factors
##
## Coef S.E. Wald Z P
## age 0.07049 0.02406 2.9304 3.386e-03
## age' -0.04335 0.02289 -1.8933 5.831e-02
## bili 0.79745 0.10922 7.3011 2.854e-13
## bili' -2.34284 0.37854 -6.1892 6.048e-10
## albumin -0.72787 0.21140 -3.4431 5.751e-04
## ascites=Yes 0.45878 0.26270 1.7464 8.074e-02
## stage=2 0.25941 0.74652 0.3475 7.282e-01
## stage=3 0.75467 0.72268 1.0443 2.964e-01
## stage=4 1.18395 0.72783 1.6267 1.038e-01
##
## Factors in Final Model
##
## [1] age bili albumin ascites stage
Compare the variables that are kept with the result for
sls=0.05. With sls=0.157 the model usually
keeps a variable or two that were removed with sls=0.05
(here ascites, which has a p-value between 0.05 and
0.157).
Why 0.157? In backward selection a variable is
removed when its p-value is larger than sls. The threshold
determines how strict the selection is:
-2 x log-likelihood + 2 x number of parameters. Removing a
variable with 1 degree of freedom lowers the AIC when the chi-square of
that variable is smaller than 2. A chi-square of 2 with 1 degree of
freedom corresponds to a p-value of 0.157. So selecting with
sls=0.157 is (for variables with 1 degree of freedom) the
same as removing variables as long as that improves the AIC. For
variables with more degrees of freedom the match is only
approximate.sls=0.05 a variable is only kept when its
chi-square is larger than 3.84 (for 1 degree of freedom). That is a much
higher bar than 2, so fewer variables are kept and the model is
smaller.A less strict threshold (0.157) keeps more variables: you are less likely to miss a variable that matters, but the model has more parameters and a larger risk of overfitting (more optimism and more shrinkage needed, see Tutorial 2). A stricter threshold (0.05) gives a smaller, simpler model, but you might remove variables that carry a bit of information. There is no single correct choice: report which threshold you used, and validate the whole procedure, including the selection, as we do in Tutorial 2. Also remember that selecting variables on p-values at all is a debatable strategy (see the note above on penalized regression).
Exercise 4. What does the cox.zph
output tell you? Which predictor (if any) might violate proportional
hazards?
survival::cox.zph(full_model)
## chisq df p
## rcs(age, 3) 0.5573 66674 1
## rcs(bili, 3) 0.4378 25006 1
## albumin 0.6281 103533 1
## sex 0.5998 203798 1
## ascites 2.2708 31013744 1
## hepato 0.0879 337198 1
## stage 4.3057 252 1
## GLOBAL 6.8567 31752713 1
Each row tests whether the effect of that predictor changes over
time. A small p-value (for example below 0.05) suggests that the
proportional hazards assumption is violated for that predictor, and the
GLOBAL row tests all predictors together. If a predictor
violates the assumption, you could add an interaction with time,
stratify on it, or use a different model. With few events these tests
have little power, so also look at the plots.
Exercise 5. Refit the full model with 4 knots for
bilirubin, rcs(bili, 4), instead of 3. Does the estimated
effect of bilirubin change? How many parameters does it add?
We fit the model in a new object, full_model_4knots, so
full_model stays as it is.
form_4knots <- S_10 ~ rcs(age, 3) + rcs(bili, 4) + albumin + sex +
ascites + hepato + stage
full_model_4knots <- Hmisc::fit.mult.impute(form_4knots,
cph,
fitargs=list(x=TRUE, y=TRUE, surv=TRUE),
xtrans=imputed_data,
data=data,
fit.reps=TRUE)
##
## Wald Statistic Information
##
## Variance Inflation Factors Due to Imputation:
##
## age age' bili bili' bili'' albumin
## 1.01 1.01 1.01 1.01 1.01 1.00
## sex=Female ascites=Yes hepato=Yes stage=2 stage=3 stage=4
## 1.00 1.14 1.09 1.08 1.07 1.08
##
## Fraction of Missing Information:
##
## age age' bili bili' bili'' albumin
## 0.01 0.01 0.01 0.01 0.01 0.00
## sex=Female ascites=Yes hepato=Yes stage=2 stage=3 stage=4
## 0.00 0.12 0.08 0.07 0.07 0.08
##
## d.f. for t-distribution for Tests of Single Coefficients:
##
## age age' bili bili' bili'' albumin
## 44243.18 20152.96 82587.68 62450.86 61604.99 26275689.48
## sex=Female ascites=Yes hepato=Yes stage=2 stage=3 stage=4
## 497354.41 267.67 571.30 741.29 881.78 669.50
##
## The following fit components were averaged over the 5 model fits:
##
## linear.predictors means stats center
# effect of bilirubin with 3 and 4 knots
plot(Predict(full_model, bili))
plot(Predict(full_model_4knots, bili))
# number of coefficients
length(coef(full_model))
## [1] 11
length(coef(full_model_4knots))
## [1] 12
With 4 knots the curve can bend in more places, but with this number of events the shape is very similar to the curve with 3 knots, while the model needs one more parameter. More knots give more flexibility, but also a higher risk of overfitting, so 3 knots is a sensible default in a data set of this size.
Exercise 6. Add edema and
copper to form_full. How does variable
importance change?
We fit the larger model in a new object,
full_model_extra, so full_model stays as it
is.
form_extra <- S_10 ~ rcs(age, 3) + rcs(bili, 3) + albumin + sex +
ascites + hepato + stage + edema + copper
full_model_extra <- Hmisc::fit.mult.impute(form_extra,
cph,
fitargs=list(x=TRUE, y=TRUE, surv=TRUE),
xtrans=imputed_data,
data=data,
fit.reps=TRUE)
##
## Wald Statistic Information
##
## Variance Inflation Factors Due to Imputation:
##
## age age'
## 1.00 1.01
## bili bili'
## 1.03 1.03
## albumin sex=Female
## 1.01 1.06
## ascites=Yes hepato=Yes
## 1.22 1.05
## stage=2 stage=3
## 1.07 1.07
## stage=4 edema=Untreated or successfully treated
## 1.07 1.02
## edema=Despite diuretics copper
## 1.01 1.79
##
## Fraction of Missing Information:
##
## age age'
## 0.00 0.01
## bili bili'
## 0.03 0.02
## albumin sex=Female
## 0.01 0.05
## ascites=Yes hepato=Yes
## 0.18 0.05
## stage=2 stage=3
## 0.07 0.07
## stage=4 edema=Untreated or successfully treated
## 0.06 0.02
## edema=Despite diuretics copper
## 0.01 0.44
##
## d.f. for t-distribution for Tests of Single Coefficients:
##
## age age'
## 257842.87 108635.55
## bili bili'
## 5498.37 6417.54
## albumin sex=Female
## 104515.53 1369.17
## ascites=Yes hepato=Yes
## 123.74 1760.96
## stage=2 stage=3
## 842.00 825.73
## stage=4 edema=Untreated or successfully treated
## 951.11 9363.50
## edema=Despite diuretics copper
## 36821.90 20.51
##
## The following fit components were averaged over the 5 model fits:
##
## linear.predictors means stats center
anova(full_model_extra)
## Wald Statistics Response: S_10
##
## Factor Chi-Square d.f. P
## age 13.18 2 0.0014
## Nonlinear 3.67 1 0.0553
## bili 66.07 2 <.0001
## Nonlinear 31.75 1 <.0001
## albumin 6.12 1 0.0133
## sex 0.00 1 0.9914
## ascites 0.82 1 0.3658
## hepato 0.11 1 0.7360
## stage 8.98 3 0.0296
## edema 13.24 2 0.0013
## copper 2.15 1 0.1429
## TOTAL NONLINEAR 32.71 2 <.0001
## TOTAL 236.75 14 <.0001
Compare the chi-square values (minus their degrees of freedom) with
the table for full_model. Copper is a known predictor in
primary biliary cirrhosis, and it is missing for about a quarter of the
patients, so its effect is estimated on imputed values. Note that
copper is right-skewed, so a transformation (for example a
log transformation) could be considered.