PREDICT - Practical 1

Getting Started with R

1 Introduction

In the next two practicals we will develop and validate prediction models. Those practicals are about prediction modelling, and they assume a certain level of understanding of R. This practical is meant to obtain and test this.

We start with the basics, but we cover them fairly quickly. The later sections deal with subsetting, factors, lists, formulas, functions, and loops. These are used constantly in Practicals 2 and 3, usually without much explanation.

How to use this practical

Reading this document is not enough. Open RStudio, create a new script, and type the code out yourself. Every exercise has a worked solution that can be unfolded, but it is worth spending a few minutes on the exercise before looking at it.

The exercises in particular are worth doing.

2 Part I: Foundations

2.1 Packages

Most functionality in R lives in a package. A package has to be installed once on your machine, and loaded again in every new R session in which you want to use it.

1install.packages("stringr")

1

Download and install the package. You do this once, ever.

There are then two ways to use a function from an installed package:

1library(stringr)
2str_to_upper("hello world")

3stringr::str_to_upper("hello world")

1

Load the whole package, making all of its functions available.

2

After loading, we can call the function directly.

3

Alternatively, package::function() calls a function without loading the package. This is more typing, but it makes it unambiguous where the function came from.

[1] "HELLO WORLD"
[1] "HELLO WORLD"

In Practicals 2 and 3 the :: form is used a lot (rms::lrm(), pROC::roc(), ggplot2::ggplot()). The reason is that when a script uses ten packages, :: makes clear which package each function comes from, and it avoids problems when two packages define a function with the same name.

Tip: you only have to install a package once, but you must load it with library() every time you (re)start R.

Tip: an Rtools installation warning on Windows can be resolved at https://cran.rstudio.com/bin/windows/Rtools/.

Practical 2 opens with a third option, pacman::p_load(), which loads several packages at once and installs any that are missing. It is therefore useful to install {pacman} now if you do not have it already installed on your computer:

install.packages("pacman")

1pacman::p_load("dplyr", "ggplot2")

1

Loads both packages, installing either of them first if it is not yet on your system.

2.2 Scripts, the console, and reproducibility

The console is useful for a quick calculation or for inspecting an object, but code typed there is lost as soon as R is closed. Anything we want to reproduce, share, or hand to a reviewer belongs in a script (.R).

A second requirement for reproducibility is set.seed(). Many of the procedures in the next practicals involve randomness, such as multiple imputation, bootstrapping, cross-validation, and splitting data. Without a seed, we obtain a different answer every time the script is run.

1set.seed(1)
2rnorm(3)

3set.seed(1)
rnorm(3)

1

Fix the state of R’s random number generator.

2

Draw three values from a standard normal distribution.

3

Re-setting the same seed reproduces exactly the same three numbers.

[1] -0.6264538  0.1836433 -0.8356286
[1] -0.6264538  0.1836433 -0.8356286

This is why the scripts in Practicals 2 and 3 all start with set.seed(1).

2.3 Loading the data

2.3.1 Traumatic brain injury (TBI) data

This dataset contains 2,159 patients from the international and US Tirilazad trials (distributed here for didactic purposes only). The primary outcome was the Glasgow Outcome (range 1 through 5) at 6 months.

The TBI dataset was sent to you together with this practical. You can also download it directly from this page:

⬇ Download TBI_data.txt

Once downloaded, place the file somewhere on your own device and load it from there:

# Specify path of TBI dataset
path <- "C:/users/avid_PREDICT_student/documents/pint/TBI.txt"

# Load TBI data
tbi <- rio::import(path)

Warning in (function (input = "", file = NULL, text = NULL, cmd = NULL, :
Detected 24 column names but the data has 25 columns (i.e. invalid file). Added
an extra default column name for the first column which is guessed to be row
names or an index. Use setnames() afterwards if this guess is not correct, or
fix the file write command that created the file to create a valid file.

Tip: make sure to use forward slashes (/) in your file paths!

The warning we receive is the result of row numbers being present in the .txt file (as also guessed by import()). Given that import() guessed correct, we can ignore this warning and treat the extra column V1 as individual identifiers.

2.3.2 Inspecting a new dataset

It is good practice to look at a dataset before doing anything with it. The four functions below are the ones we use most often for this.

1dim(tbi)
2names(tbi)
3head(tbi, n = 3)
4str(tbi[, 1:8])

1

Number of rows and columns.

2

The column names.

3

The first three rows.

4

The structure: for each column, its type and first few values. str() is probably the most useful function in R for checking whether an object contains what we think it does.

[1] 2159   25
 [1] "V1"       "trial"    "d.gos"    "d.mort"   "d.unfav"  "cause"   
 [7] "age"      "d.motor"  "d.pupil"  "pupil.i"  "hypoxia"  "hypotens"
[13] "ctclass"  "tsah"     "edh"      "cisterns" "shift"    "d.sysbpt"
[19] "glucose"  "glucoset" "ph"       "sodium"   "sodiumt"  "hb"      
[25] "hbt"     
  V1                   trial         d.gos d.mort d.unfav         cause age
1  1 Tirilazad International good recovery      0       0     Motorbike  14
2  2 Tirilazad International good recovery      0       0     Motorbike  14
3  3 Tirilazad International good recovery      0       0 domestic/fall  14
  d.motor       d.pupil       pupil.i hypoxia hypotens ctclass tsah edh
1       5 both reactive both reactive       0        0       2    0   0
2       4 both reactive both reactive       0        0       2    0   0
3       4 both reactive both reactive       1        0       4    1   0
  cisterns shift d.sysbpt  glucose glucoset    ph sodium sodiumt   hb  hbt
1        1     0   119.09 7.700000 7.700000 7.350    143     143 15.0 15.0
2        1     0   130.71 7.400000 7.400000 7.330    143     143 15.0 15.0
3       NA    NA   136.26 7.555556 7.555556 7.489    141     141 12.8 12.8
'data.frame':   2159 obs. of  8 variables:
 $ V1     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ trial  : chr  "Tirilazad International" "Tirilazad International" "Tirilazad International" "Tirilazad International" ...
 $ d.gos  : chr  "good recovery" "good recovery" "good recovery" "moderate disability" ...
 $ d.mort : int  0 0 0 0 0 0 0 0 0 0 ...
 $ d.unfav: int  0 0 0 0 0 0 0 0 0 0 ...
 $ cause  : chr  "Motorbike" "Motorbike" "domestic/fall" "Motorbike" ...
 $ age    : int  14 14 14 14 14 14 14 14 14 14 ...
 $ d.motor: int  5 4 4 4 5 3 5 5 4 5 ...

2.4 Vectors, the basic building block

The c() function combines values into a vector. Most objects in R are built out of vectors.

1a <- c(1, 2, 3)
2b <- 1:3
3d <- seq(0, 10, by = 2.5)
4e <- rep(1, 5)

5length(d)
6a * 2
7a + b

1

Combine three values into a vector and store it in a.

2

1:3 is shorthand for a sequence from 1 to 3.

3

seq() builds a sequence with a chosen step size.

4

rep() repeats a value; here, five ones.

5

Counts the number of elements the vector has.

6

Arithmetic is vectorised, meaning that the operation is applied to every element at once, without a loop.

7

Two vectors of the same length are combined element by element.

[1] 5
[1] 2 4 6
[1] 2 4 6

Vectorisation

R is built around operations on whole vectors. For example, tbi$age * 12 converts 2,159 ages into months in one step. If we find ourselves writing a loop to do something to each element of a vector, there is usually a vectorised alternative that is both shorter and faster.

2.5 Types and coercion

A vector holds exactly one type of data. The three types we encounter in this course are numeric, character, and logical.

1class(tbi$age)
2class(tbi$cause)
3class(tbi$age > 50)

1

age holds numbers.

2

cause holds text such as "Assault" or "Motorbike".

3

A comparison returns TRUE/FALSE values: a logical vector.

[1] "integer"
[1] "character"
[1] "logical"

Because a vector can only hold one type, R silently coerces everything to the most flexible type present when types are mixed:

1c(1, 2, 3)
2c(1, 2, "3")
3c(1, 2, TRUE)

1

All numeric.

2

One character value forces the whole vector to become character. Note the quotes around every element.

3

TRUE becomes 1; logicals sit “below” numbers in the hierarchy.

[1] 1 2 3
[1] "1" "2" "3"
[1] 1 2 1

The last rule is more useful than it may look:

1sum(tbi$age > 50)
2mean(tbi$age > 50)
3as.numeric(tbi$cause == "Assault") |> head()

1

TRUE counts as 1 and FALSE as 0, so summing a logical vector counts how many are TRUE.

2

Averaging a logical vector gives the proportion that are TRUE.

3

Explicit conversion to 0/1, which is often what a model wants.

[1] 299
[1] 0.13849
[1] 0 0 0 0 0 0

We can convert types on purpose with the as.*() family, but only when the conversion makes sense:

1as.numeric(tbi$cause) |> head()

1

R cannot read "Assault" as a number, so it returns NA and warns you. Whenever you see “NAs introduced by coercion”, you have tried to turn text into numbers.

[1] NA NA NA NA NA NA

Dates are a type of their own and can behave in unexpected ways, so it is worth checking the class of a date variable before doing arithmetic with it. Adding a year to a date, for example, requires a package such as {lubridate}:

Sys.Date() + lubridate::years(1)

2.6 Missing and undefined values

R distinguishes three things that all look like “nothing”:

1sum(is.na(tbi$d.pupil))
20 / 0
3tbi$region

1

NA is a missing value: the patient exists, the value is unknown. is.na() returns a logical vector, and summing it counts the missings.

2

NaN (“not a number”) is the result of an undefined mathematical operation.

3

NULL is the absence of an object. tbi has no column region, so R returns NULL.

[1] 123
[1] NaN
NULL

Missing values propagate through calculations, which is a common source of confusion:

x <- c(1, 2, NA, 4)

1mean(x)
2mean(x, na.rm = TRUE)
3x == NA
4is.na(x)

1

The mean of something partly unknown is unknown: NA.

2

na.rm = TRUE tells the function to drop missings first. Many summary functions have this argument.

3

This does not work. The question “is this unknown value equal to an unknown value?” is itself unknown.

4

This is the correct way to test for missingness. Always use is.na().

[1] NA
[1] 2.333333
[1] NA NA NA NA
[1] FALSE FALSE  TRUE FALSE

Exercise 1

Using the tbi data:

1. How many patients are older than 65?
2. What percentage of patients had an unfavourable outcome (d.unfav)? Do this in one line, without using table().
3. What is the mean glucose (glucose) at admission? Check first whether the variable contains missings.
4. Create a logical vector that is TRUE for patients who are both older than 65 and had an unfavourable outcome, and count them.

Solution 1

# a. TRUE counts as 1, so summing the logical vector counts patients
sum(tbi$age > 65)

[1] 20

# b. the mean of a logical vector is the proportion; times 100 for a percentage
mean(tbi$d.unfav) * 100

[1] 39.4164

# c. first check for missings, then compute
sum(is.na(tbi$glucose))

[1] 64

mean(tbi$glucose, na.rm = TRUE)

[1] 8.907855

# d. & combines two logical vectors element by element
old_and_unfav <- tbi$age > 65 & tbi$d.unfav == 1
sum(old_and_unfav)

[1] 16

3 Part II: Data structures

Practicals 2 and 3 spend much of their time getting things out of objects: columns out of data frames, coefficients out of model fits, and levels out of categorical variables. This part covers what is needed for that, namely subsetting, factors, lists, and formulas.

3.1 Matrices and data frames

A matrix is rectangular and holds a single type throughout, which makes it suitable for mathematical operations.

1mat <- matrix(1:6, nrow = 2, byrow = TRUE)
mat
2dim(mat)
3mat * mat
4mat %*% t(mat)

1

A 2-by-3 matrix, filled by row.

2

Two rows, three columns.

3

Element-wise multiplication.

4

True matrix multiplication with the transpose, using %*%.

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[1] 2 3
     [,1] [,2] [,3]
[1,]    1    4    9
[2,]   16   25   36
     [,1] [,2]
[1,]   14   32
[2,]   32   77

Matrix algebra reappears in Practical 3, where the model matrix x is multiplied by a coefficient vector to recompute linear predictors by hand.

A data frame is also rectangular, but each column may have its own type. This is the standard structure for patient data, and it is what tbi is.

1patients <- data.frame(
  id   = 1:3,
  name = c("Anna", "Ben", "Chen"),
  age  = c(28, 34, 51)
)
patients
class(patients)
2sapply(patients, class)

1

A data frame with one integer, one character, and one numeric column.

2

Ask for the class of each column at once. (More on sapply() later.)

  id name age
1  1 Anna  28
2  2  Ben  34
3  3 Chen  51
[1] "data.frame"
         id        name         age 
  "integer" "character"   "numeric" 

3.2 Subsetting: [, [[, and $

R has three subsetting operators, and a fair amount of confusion comes from mixing them up.

3.2.1 Subsetting a data frame with [row, column]

1tbi[1, ]
2tbi[1:3, c("age", "cause")]
3tbi[, "age"] |> head()
4tbi[tbi$age > 70, c("age", "d.unfav")] |> head()
5tbi[-(1:2150), c("age", "cause")]

1

Row 1, all columns. The comma is essential: tbi[1, ] is a row, tbi[1] is a column.

2

Rows 1 to 3, and only the columns age and cause.

3

Leaving the row slot empty takes all rows. This returns the column as a plain vector.

4

Rows can also be selected with a logical vector, which is what filtering does internally.

5

A negative index means “everything except”. Here we drop the first 2,150 rows.

  V1                   trial         d.gos d.mort d.unfav     cause age d.motor
1  1 Tirilazad International good recovery      0       0 Motorbike  14       5
        d.pupil       pupil.i hypoxia hypotens ctclass tsah edh cisterns shift
1 both reactive both reactive       0        0       2    0   0        1     0
  d.sysbpt glucose glucoset   ph sodium sodiumt hb hbt
1   119.09     7.7      7.7 7.35    143     143 15  15
  age         cause
1  14     Motorbike
2  14     Motorbike
3  14 domestic/fall
[1] 14 14 14 14 14 14
     age d.unfav
1113  71       1
1114  73       1
1115  74       0
1116  75       0
1117  76       1
1118  79       1
     age         cause
2151  65 domestic/fall
2152  65       Assault
2153  68         other
2154  68 domestic/fall
2155  68 domestic/fall
2156  68     Motorbike
2157  69         other
2158  70 domestic/fall
2159  77 domestic/fall

3.2.2 [ versus [[

The rule is that [ keeps the container, whereas [[ takes a single element out of it.

1class(tbi["age"])
2class(tbi[["age"]])
3class(tbi$age)

1

A data frame with one column.

2

The column itself, as a numeric vector.

3

$ does the same as [[ here, with less typing.

[1] "data.frame"
[1] "integer"
[1] "integer"

Practicals 2 and 3 nearly always write tbi[["age"]] rather than tbi$age. There are two reasons for this.

1tbi$ag |> head()

1

$ performs partial matching, so ag is silently resolved to age. This is convenient when typing in the console, but in a script it means that a typo returns the wrong column instead of an error.

[1] 14 14 14 14 14 14

my_var <- "age"
1tbi$my_var
2tbi[[my_var]] |> head()

1

$ takes the name literally, looks for a column called my_var, does not find one, and returns NULL.

2

[[ evaluates the variable first, so it does find the age column. This is what makes loops over variable names possible, which we need in Practical 3.

NULL
[1] 14 14 14 14 14 14

Which one to use

$ is fine when typing in the console. In a script, [["name"]] is safer, and it is the only option when the column name is stored in a variable.

3.2.3 Adding and removing columns

1tbi$assault <- as.numeric(tbi$cause == "Assault")
2tbi[["elderly"]] <- as.numeric(tbi$age > 65)
tbi[, c("assault", "elderly")] |> head()

3tbi$assault <- NULL
4"assault" %in% names(tbi)

1

Assigning to a column name that does not exist creates it.

2

The [[ form does exactly the same and is preferred in scripts.

3

Assigning NULL deletes a column.

4

%in% tests membership and returns TRUE/FALSE. You will see "glmnet" %in% class(fit) in Practical 3.

  assault elderly
1       0       0
2       0       0
3       0       0
4       0       0
5       0       0
6       0       0
[1] FALSE

Exercise 2

1. Select the columns age, d.motor, and d.unfav for the first five patients.
2. Extract the age column in three different ways, and confirm with class() that two of them give a vector and one gives a data frame.
3. Create a subset containing only patients from the "Tirilazad US" trial. How many rows does it have?
4. Store the string "d.motor" in a variable called v. Now compute the mean of that column using v, without typing d.motor again.
5. Add a column young that is 1 for patients under 30 and 0 otherwise, then remove it again.

Solution 2

# a.
tbi[1:5, c("age", "d.motor", "d.unfav")]

  age d.motor d.unfav
1  14       5       0
2  14       4       0
3  14       4       0
4  14       4       0
5  14       5       0

# b.
class(tbi$age)        # vector

[1] "integer"

class(tbi[["age"]])   # vector

[1] "integer"

class(tbi["age"])     # data frame

[1] "data.frame"

# c.
us <- tbi[tbi$trial == "Tirilazad US", ]
nrow(us)

[1] 1041

# d. only [[ ]] can look up a name stored in a variable
v <- "d.motor"
mean(tbi[[v]])

[1] 3.9912

# e.
tbi[["young"]] <- as.numeric(tbi$age < 30)
table(tbi[["young"]])

   0    1 
1103 1056 

tbi[["young"]] <- NULL

3.3 Factors: how R stores categorical variables

A factor is the type R uses for categorical variables. It looks like text, but it is stored as integer codes together with a set of levels. This is not only a matter of presentation, because it changes the model that is fitted.

1cause_f <- factor(tbi$cause)
class(cause_f)
2levels(cause_f)
3table(cause_f)
4as.integer(cause_f) |> head()

1

Convert the character column into a factor.

2

The categories, in alphabetical order by default.

3

How many patients in each category.

4

Underneath, a factor really is a vector of integers pointing at the levels.

[1] "factor"
[1] "Assault"               "domestic/fall"         "Motorbike"            
[4] "other"                 "Road traffic accident"
cause_f
              Assault         domestic/fall             Motorbike 
                  134                   370                   420 
                other Road traffic accident 
                  387                   848 
[1] 3 3 2 3 3 3

3.3.1 Why it matters for modelling

Take d.motor, the admission motor score, which is imported as a number from 1 to 6.

1fit_num <- glm(d.unfav ~ d.motor, data = tbi, family = binomial)
coef(fit_num)

2fit_fac <- glm(d.unfav ~ factor(d.motor), data = tbi, family = binomial)
coef(fit_fac)

1

Treated as numeric, d.motor gets one coefficient. The model then assumes that the step from 1 to 2 has the same effect as the step from 5 to 6.

2

Treated as a factor, it gets five coefficients, one for each category other than the first, and these are free to differ. Category 1 has become the reference category.

(Intercept)     d.motor 
  2.2739937  -0.6886543 
     (Intercept) factor(d.motor)2 factor(d.motor)3 factor(d.motor)4 
       0.5877867        0.3410828       -0.3646431       -1.0386534 
factor(d.motor)5 factor(d.motor)6 
      -1.8082888       -2.2275299 

Always declare your categorical variables as factors

If a categorical variable is coded with numbers, R has no way of knowing that it is categorical. It treats the variable as continuous and fits a single coefficient, which is usually not meaningful. Practical 2 performs this conversion before modelling:

tbi[["d.motor"]] <- as.factor(tbi[["d.motor"]])

3.3.2 Controlling levels, labels, and the reference category

The first level is the reference category, so its order matters for interpretation.

pupil_f <- factor(
  tbi$d.pupil,
1  levels = c("no reactive pupils", "one reactive", "both reactive"),
2  labels = c("none", "one", "both")
)
3table(pupil_f, useNA = "ifany")
4levels(pupil_f)[1]

5pupil_ref <- relevel(pupil_f, ref = "both")
levels(pupil_ref)

1

levels states which values exist and in what order. We choose the order ourselves, and it is only alphabetical if we make it so.

2

labels renames them, in the same order. This is how Practical 2 turns the pupil text into 0, 1, 2.

3

useNA = "ifany" also shows the missing values, which table() hides by default. It is worth checking this at least once for every categorical variable.

4

The first level is the reference category against which the coefficients are compared.

5

relevel() makes a different category the reference. This changes the coefficients, but not the fit of the model itself.

pupil_f
none  one both <NA> 
 327  279 1430  123 
[1] "none"
[1] "both" "none" "one" 

The classic factor trap

If the labels of a factor are numbers, converting it with as.numeric() returns the internal codes rather than the numbers we see.

f <- factor(c(10, 20, 30))
1as.numeric(f)
2as.numeric(as.character(f))

1

Wrong: you get the level codes 1, 2, 3.

2

Right: go via as.character() first.

[1] 1 2 3
[1] 10 20 30

Exercise 3

1. Turn cause into a factor whose reference category is "Road traffic accident", with the remaining categories in any order you like.
2. Fit glm(d.unfav ~ cause_factor, ...) with family = binomial and count how many coefficients it has. Explain the number.
3. Refit with "Assault" as the reference category using relevel(). Which numbers change, and which do not? (Compare the model deviances with deviance().)

Solution 3

# a.
tbi$cause_f <- factor(
  tbi$cause,
  levels = c("Road traffic accident", "Motorbike",
             "domestic/fall", "Assault", "other")
)
levels(tbi$cause_f)

[1] "Road traffic accident" "Motorbike"             "domestic/fall"        
[4] "Assault"               "other"                

# b. five categories -> one intercept + four contrasts = 5 coefficients
fit_a <- glm(d.unfav ~ cause_f, data = tbi, family = binomial)
coef(fit_a)

         (Intercept)     cause_fMotorbike cause_fdomestic/fall 
          -0.5108256           -0.1503473            0.3483067 
      cause_fAssault         cause_fother 
          -0.1049349            0.2981412 

length(coef(fit_a))

[1] 5

# c. the coefficients change because they are compared to a different
#    reference, but the model itself is identical: same deviance.
tbi$cause_ref <- relevel(tbi$cause_f, ref = "Assault")
fit_c <- glm(d.unfav ~ cause_ref, data = tbi, family = binomial)
coef(fit_c)

                   (Intercept) cause_refRoad traffic accident 
                   -0.61576052                     0.10493489 
            cause_refMotorbike         cause_refdomestic/fall 
                   -0.04541236                     0.45324159 
                cause_refother 
                    0.40307610 

c(deviance(fit_a), deviance(fit_c))

[1] 2877.023 2877.023

3.4 Lists, and why a model fit is one

A list is a container that can hold anything, including objects of different types and lengths. It is the most flexible structure in R, and most complicated R objects are built out of lists.

1info <- list(
  study     = "Tirilazad",
  n         = nrow(tbi),
  predictors = c("age", "d.motor", "d.pupil")
)

2names(info)
3info[["n"]]
4info$predictors[2]
5class(info["n"])
6class(info[["n"]])
7str(info)

1

One character value, one number, and a character vector of length three, in a single object.

2

The names of the elements.

3

Pull out one element with [[.

4

Elements can themselves be vectors, so we can keep indexing.

5

[ returns a list containing that element.

6

[[ returns the element itself. This is the same rule as for data frames, which are in fact a special kind of list.

7

str() shows the shape of the whole thing.

[1] "study"      "n"          "predictors"
[1] 2159
[1] "d.motor"
[1] "list"
[1] "integer"
List of 3
 $ study     : chr "Tirilazad"
 $ n         : int 2159
 $ predictors: chr [1:3] "age" "d.motor" "d.pupil"

3.4.1 The model fit object

When we fit a model, what we get back is a list.

1fit <- glm(d.unfav ~ age + factor(d.motor), data = tbi, family = binomial)

2class(fit)
3names(fit)

1

Fit a logistic regression model.

2

Its class is glm (and lm), which is what tells summary() and predict() how to behave.

3

But underneath it is just a named list, and these are its elements.

[1] "glm" "lm" 
 [1] "coefficients"      "residuals"         "fitted.values"    
 [4] "effects"           "R"                 "rank"             
 [7] "qr"                "family"            "linear.predictors"
[10] "deviance"          "aic"               "null.deviance"    
[13] "iter"              "weights"           "prior.weights"    
[16] "df.residual"       "df.null"           "y"                
[19] "converged"         "boundary"          "model"            
[22] "call"              "formula"           "terms"            
[25] "data"              "offset"            "control"          
[28] "method"            "contrasts"         "xlevels"          

Most of what Practical 2 does with a model comes down to taking elements out of this list:

1fit[["coefficients"]]
2fit[["coefficients"]][["age"]]
3coef(fit)[["(Intercept)"]]
4head(fit[["fitted.values"]])
5head(fit[["linear.predictors"]])
6fit[["df.residual"]]

1

All estimated coefficients, as a named vector.

2

Named vectors can be indexed by name with [[ too. This gives the log odds ratio for age.

3

Most model objects also have extractor functions that do the same job in a safer way. The common ones are coef(), residuals(), fitted(), and predict().

4

The predicted probabilities for every patient. Practical 2 assigns exactly this to a new column.

5

The linear predictors (the log odds), which Practical 2 uses to compute the calibration intercept and slope.

6

Residual degrees of freedom.

     (Intercept)              age factor(d.motor)2 factor(d.motor)3 
     -0.74859525       0.03749433       0.50441704      -0.22763980 
factor(d.motor)4 factor(d.motor)5 factor(d.motor)6 
     -0.99526136      -1.77747771      -2.26360472 
[1] 0.03749433
[1] -0.7485952
        1         2         3         4         5         6 
0.1190820 0.2281238 0.2281238 0.2281238 0.1190820 0.3890483 
         1          2          3          4          5          6 
-2.0011523 -1.2189359 -1.2189359 -1.2189359 -2.0011523 -0.4513144 
[1] 2152

When you do not know what is in an object

names(), str(), and class() will tell us. If str(fit) prints too much, str(fit, max.level = 1) gives a shorter overview. This is usually the quickest way to find out what an unfamiliar model object contains.

Because the model’s class determines how functions behave, code that must handle several model types tests it explicitly:

1class(fit)
2"glm" %in% class(fit)
3inherits(fit, "glm")

1

An object can have more than one class.

2

%in% checks whether "glm" is one of them, and returns a single TRUE/FALSE.

3

inherits() is the tidier way of writing the same test.

[1] "glm" "lm" 
[1] TRUE
[1] TRUE

Exercise 4

Fit the model d.unfav ~ age + factor(d.pupil) on the tbi data with family = binomial.

1. How many elements does the resulting object contain?
2. Extract the coefficient for age in two different ways.
3. Extract the predicted probabilities and compute their mean. Compare it with the observed proportion of unfavourable outcomes among the patients used in the fit. (Hint: fit[["y"]] holds the outcomes that were actually used.)
4. Why are fewer patients used than the 2,159 in the dataset?

Solution 4

fit_ex <- glm(d.unfav ~ age + factor(d.pupil),
              data = tbi, family = binomial)

# a.
length(fit_ex)

[1] 31

# b.
fit_ex[["coefficients"]][["age"]]

[1] 0.03433832

coef(fit_ex)[["age"]]

[1] 0.03433832

# c. for a correctly fitted logistic model these agree almost exactly
mean(fit_ex[["fitted.values"]])

[1] 0.391945

mean(fit_ex[["y"]])

[1] 0.391945

# d. d.pupil has missing values, and glm() drops those patients by default
sum(is.na(tbi$d.pupil))

[1] 123

nrow(tbi) - length(fit_ex[["y"]])

[1] 123

Part (d) is important, because R performs a complete-case analysis here. The summary of the model does mention this, but only in a single line at the bottom, which is easy to miss. Practical 2 deals with it using imputation.

3.5 Formulas

Every model in Practicals 2 and 3 is specified with a formula of the form outcome ~ predictors. A formula is an object in its own right, and not just part of the function call.

1f <- d.unfav ~ age + d.motor + cause
2class(f)
3all.vars(f)

1

Read the tilde as “is modelled as a function of”.

2

Formulas have their own class; they can be stored, passed around, and modified.

3

Which variables appear in it.

[1] "formula"
[1] "d.unfav" "age"     "d.motor" "cause"  

Because a formula is an object, we can define a model in one place and use it in another:

1fit_f <- glm(f, data = tbi, family = binomial)
2length(coef(fit_f))

3f2 <- update(f, . ~ . + hypoxia)
f2

1

The formula is passed to glm() exactly like any other argument.

2

cause is still a character variable here. glm() converts character predictors into factors automatically, which is why we obtain more than three coefficients. It is better not to rely on this and to declare the factors ourselves.

3

update() modifies a formula. The dots mean “keep the existing outcome” and “keep the existing predictors”, and we add one more.

[1] 7
d.unfav ~ age + d.motor + cause + hypoxia

Some formula notation that we will come across:

Formula notation you will encounter

Notation	Meaning
y ~ x1 + x2	Two predictors, additive
y ~ 1	Intercept only, no predictors
y ~ .	All remaining columns as predictors
y ~ x1 * x2	Both main effects and their interaction
y ~ offset(lp)	Include lp with its coefficient fixed at 1
rms::rcs(age, 4)	Age as a restricted cubic spline with 4 knots
Surv(time, status) ~ x	A survival outcome, for Cox models

The offset() term is worth a separate note, because it is not obvious what it does when it first appears in Practical 2. Fitting y ~ offset(lp) estimates only an intercept, while the coefficient of lp is fixed at 1. That is exactly the definition of the calibration intercept. Below you find a preview on how to compute this, but we will go into more detail in Practical 2.

1cal_data <- data.frame(
  d.unfav = fit[["y"]],
  lp      = fit[["linear.predictors"]]
)

2cal_int <- coef(glm(d.unfav ~ offset(lp),
                    data = cal_data, family = binomial))[[1]]
3cal_slope <- coef(glm(d.unfav ~ lp,
                      data = cal_data, family = binomial))[["lp"]]

4c(intercept = cal_int, slope = cal_slope)

1

Collect the outcome and the linear predictor in a small data frame of their own. Taking the outcome from fit[["y"]] rather than from tbi ensures that the two have the same length, also when the model dropped patients with missing values.

2

The calibration intercept: the only free parameter is the intercept, because offset() fixes the coefficient of lp at 1.

3

The calibration slope: the coefficient of the linear predictor.

4

Evaluated on the same data used for development, these are 0 and 1 by construction. In Practical 2 you will see why that is the expected answer, and what it means when new data gives something else.

   intercept        slope 
2.506454e-11 1.000000e+00 

4 Part III: Working with data

4.1 Filtering and modifying

Rows can be filtered with base R’s subset(), or with dplyr::filter(), which most people find more readable in a chain.

library(dplyr)

1subset(tbi, age > 50) |> nrow()
2filter(tbi, age > 50) |> nrow()
3tbi |> filter(age > 50) |> nrow()
4tbi |> filter(age > 50, d.unfav == 1) |> nrow()

1

Base R.

2

The {dplyr} equivalent.

3

The native pipe |> passes the left-hand side as the first argument of the right-hand function. Practical 2 also uses %>%, the {magrittr} pipe, which behaves the same way here.

4

Several conditions separated by commas are combined with “and”.

[1] 299
[1] 299
[1] 299
[1] 167

mutate() creates or changes columns, and ifelse() is the function we normally use for recoding.

tbi <- tbi |>
  mutate(
1    age_months = age * 12,
2    elderly    = ifelse(age >= 65, 1, 0),
3    severity   = case_when(
      d.motor <= 2 ~ "severe",
      d.motor <= 4 ~ "moderate",
      TRUE         ~ "mild"
    )
  )

table(tbi$severity)

1

A straightforward transformation.

2

ifelse(test, yes, no) is vectorised, so it evaluates the condition for every row and picks the corresponding value.

3

case_when() handles more than two categories. Conditions are checked in order, and TRUE ~ ... is the catch-all “everything else”.

    mild moderate   severe 
     870      996      293 

Warning

ifelse() propagates NA, so if the test is NA for a patient, the result is NA as well. That is usually the desired behaviour, but it is worth checking rather than assuming.

Practical 2 uses exactly this pattern to administratively censor survival times at two years:

lung <- lung |>
    mutate(status = ifelse(time > 365.25 * 2, 1, status),
           time   = ifelse(time > 365.25 * 2, 365.25 * 2, time))

4.2 Summarising

The functions below are the ones used in the next practicals to describe the data before modelling.

1table(tbi$d.unfav)
2proportions(table(tbi$d.unfav)) * 100
3table(tbi$d.pupil, tbi$d.unfav)
4summary(tbi[, c("age", "d.motor", "glucose")])
5colSums(is.na(tbi[, c("d.unfav", "age", "d.pupil", "cause")]))

1

Counts per category.

2

The same as percentages. Practical 2 uses this to judge the risk of overfitting from the number of events.

3

A cross-tabulation of two variables.

4

A quick numerical summary of several columns at once, including a count of NAs.

5

The number of missing values per column, which we should check before fitting anything.

   0    1 
1308  851 

      0       1 
60.5836 39.4164 
                    
                        0    1
  both reactive      1008  422
  no reactive pupils   98  229
  one reactive        132  147
      age           d.motor         glucose       
 Min.   :14.00   Min.   :1.000   Min.   : 0.5556  
 1st Qu.:22.00   1st Qu.:3.000   1st Qu.: 6.7000  
 Median :30.00   Median :4.000   Median : 8.2000  
 Mean   :33.21   Mean   :3.991   Mean   : 8.9078  
 3rd Qu.:43.00   3rd Qu.:5.000   3rd Qu.:10.4000  
 Max.   :79.00   Max.   :6.000   Max.   :41.4000  
                                 NA's   :64       
d.unfav     age d.pupil   cause 
      0       0     123       0 

To compute something within groups, use by() or tapply():

1tapply(tbi$age, tbi$trial, mean)
2by(tbi$d.unfav, tbi$trial, \(x) round(proportions(table(x)) * 100, 1))

1

Mean age per trial.

2

The event percentage per trial. Here \(x) is an anonymous function, that is, a function defined on the spot without a name. We return to this below. This line also appears in Practical 3.

Tirilazad International            Tirilazad US 
               33.61270                32.78002 
tbi$trial: Tirilazad International
x
   0    1 
59.2 40.8 
------------------------------------------------------------ 
tbi$trial: Tirilazad US
x
   0    1 
62.1 37.9 

Exercise 5

1. Make a table of cause showing percentages rather than counts, rounded to one decimal.
2. Cross-tabulate d.pupil against d.unfav, including the missing values.
3. Add a column age_cat with categories "<40", "40-64", and "65+", then tabulate it against d.unfav.
4. Compute the mean age of patients with and without an unfavourable outcome, in one line.
5. Which columns of tbi contain missing values? (Hint: apply is.na() and colSums() to the whole data frame, then keep the non-zero entries.)

Solution 5

# a.
round(proportions(table(tbi$cause)) * 100, 1)

              Assault         domestic/fall             Motorbike 
                  6.2                  17.1                  19.5 
                other Road traffic accident 
                 17.9                  39.3 

# b.
table(tbi$d.pupil, tbi$d.unfav, useNA = "ifany")

                    
                        0    1
  both reactive      1008  422
  no reactive pupils   98  229
  one reactive        132  147
  <NA>                 70   53

# c.
tbi <- tbi |>
  mutate(age_cat = case_when(age < 40  ~ "<40",
                             age < 65  ~ "40-64",
                             TRUE      ~ "65+"))
table(tbi$age_cat, tbi$d.unfav)

       
           0    1
  <40   1015  498
  40-64  283  332
  65+     10   21

# d.
tapply(tbi$age, tbi$d.unfav, mean)

       0        1 
30.92966 36.71798 

# e. count missings per column, then keep only those above zero
miss <- colSums(is.na(tbi))
miss[miss > 0]

 d.pupil  hypoxia hypotens  ctclass     tsah      edh cisterns    shift 
     123      249       59       25       96       38      257      252 
 glucose glucoset       ph   sodium  sodiumt       hb      hbt 
      64       64      457       62       62       24       24 

5 Part IV: Programming

5.1 Writing functions

A function bundles a set of operations so that they can be reused. We have used plenty of them already, and now we will write some ourselves.

1say_hello <- function() {
  print("Hello world!")
}
say_hello()

2odds <- function(p) {
3  return(p / (1 - p))
}
odds(0.25)
4odds(c(0.1, 0.5, 0.9))

1

A function with no arguments.

2

A function with one argument.

3

return() states explicitly what comes out. Without it, R returns the last expression that was evaluated, but writing it out is clearer.

4

Because the body of the function is vectorised, the function itself is vectorised as well, without any extra work.

[1] "Hello world!"
[1] 0.3333333
[1] 0.1111111 1.0000000 9.0000000

5.1.1 Default arguments

Arguments can be given default values, which makes a function easier to call while keeping it flexible.

1describe <- function(x, digits = 1, na.rm = TRUE) {
  c(mean = round(mean(x, na.rm = na.rm), digits),
    sd   = round(sd(x, na.rm = na.rm), digits))
}

2describe(tbi$age)
3describe(tbi$age, digits = 3)

1

x is required; digits and na.rm have defaults.

2

Called with the defaults.

3

Overriding one of the defaults. Naming arguments explicitly is good practice, because the code then still works if the order of the arguments changes.

mean   sd 
33.2 13.6 
  mean     sd 
33.211 13.578 

5.1.2 Conditionals inside functions

if and else let a function change its behaviour depending on its input. This is how the helper functions in Practical 3 deal with two different kinds of model object.

1risk_group <- function(p, threshold = 0.5) {
2  if (any(is.na(p))) {
    warning("Missing predictions were found.")
  }
3  if (length(threshold) != 1) {
    stop("`threshold` must be a single number.")
  }
4  ifelse(p >= threshold, "high", "low")
}

risk_group(c(0.2, 0.8, 0.55))

1

One required argument and one with a default.

2

warning() prints a message but lets the function continue.

3

stop() aborts with an error. Checking the input at the start of a function usually saves time later on.

4

Note the difference between the two. if handles a single TRUE or FALSE and decides which code is run, whereas ifelse() works element by element on a vector and returns a vector.

[1] "low"  "high" "high"

5.1.3 Returning several things at once

A function can return only one object, but that object can be a list.

1cal_metrics <- function(fit, newdata = NULL) {

2  if (is.null(newdata)) {
    lp <- fit[["linear.predictors"]]
    y  <- fit[["y"]]
  } else {
3    lp <- predict(fit, newdata = newdata, type = "link")
    y  <- newdata[["d.unfav"]]
  }

4  cal_int   <- coef(glm(y ~ offset(lp), family = binomial))[[1]]
  cal_slope <- coef(glm(y ~ lp, family = binomial))[[2]]

5  return(list(intercept = cal_int, slope = cal_slope))
}

6out <- cal_metrics(fit)
out
7out[["slope"]]

1

newdata = NULL is the idiom for an optional argument: “use the development data unless I give you something else”.

2

is.null() tests whether the argument was supplied, and the two branches take the values from different places.

3

predict() with newdata applies the fitted model to patients that were not used to fit it, which is what external validation comes down to.

4

The calibration intercept and slope, exactly as in the formula section.

5

Wrapping both results in a named list lets us return them together.

6

Called without newdata, so it evaluates the model on its own development data.

7

Elements come back out with [[, as always.

$intercept
[1] 2.506454e-11

$slope
[1] 1

[1] 1

The function calmetrics() in Practical 3 has the same structure, with the same NULL defaults, the same conditional, and a list as output.

5.1.4 Anonymous functions

A function that is used only once does not need a name. The \(x) syntax is shorthand for function(x).

1sapply(c(1, 4, 9), \(x) sqrt(x))
2sapply(c(1, 4, 9), function(x) sqrt(x))

1

The modern, compact form. Practical 3 uses this inside sapply() and by().

2

Identical, written out in full.

[1] 1 2 3
[1] 1 2 3

Exercise 6

1. Write a function logit(p) returning \(\log(p / (1-p))\), and invlogit(x) returning \(1 / (1 + e^{-x})\). Check that invlogit(logit(0.3)) gives 0.3 back. (R also provides these as qlogis() and plogis().)
2. Extend logit() so that it stops with an informative error if any value is outside the interval (0, 1).
3. Write describe_model(fit) that returns a list with the number of coefficients, the number of patients used, and the mean predicted probability.
4. Run it on a model of your choice and extract just the mean predicted probability.

Solution 6

# a.
logit    <- function(p) log(p / (1 - p))
invlogit <- function(x) 1 / (1 + exp(-x))
invlogit(logit(0.3))

[1] 0.3

# b. check the input before doing anything with it
logit <- function(p) {
  if (any(p <= 0 | p >= 1, na.rm = TRUE)) {
    stop("All values of `p` must lie strictly between 0 and 1.")
  }
  log(p / (1 - p))
}
logit(0.3)

[1] -0.8472979

# c.
describe_model <- function(fit) {
  return(list(
    n_coef  = length(coef(fit)),
    n_used  = length(fit[["y"]]),
    mean_p  = mean(fit[["fitted.values"]])
  ))
}

# d.
info <- describe_model(fit)
str(info)

List of 3
 $ n_coef: int 7
 $ n_used: int 2159
 $ mean_p: num 0.394

info[["mean_p"]]

[1] 0.394164

5.2 Loops and the apply family

5.2.1 For-loops

A for-loop repeats a block of code once for each element of a sequence.

1for (i in 1:3) {
  cat("iteration", i, "\n")
}

2for (v in c("age", "glucose")) {
  cat(v, "has mean", round(mean(tbi[[v]], na.rm = TRUE), 2), "\n")
}

1

i takes the values 1, 2, 3 in turn.

2

The loop variable does not have to be a number. Here it holds column names, which we then look up with [[. This works with [[ but not with $, which is why we spent some time on the difference earlier.

iteration 1 
iteration 2 
iteration 3 
age has mean 33.21 
glucose has mean 8.91 

5.2.2 Collecting results

A loop that only prints is not very useful, since we normally want to keep the results. There are two common patterns for doing this.

vars <- c("age", "d.motor", "glucose")

1means <- numeric(length(vars))
names(means) <- vars
2for (i in seq_along(vars)) {
  means[i] <- mean(tbi[[vars[i]]], na.rm = TRUE)
}
round(means, 2)

1

Pre-allocate, that is, first create an empty container of the right size. This is the more efficient of the two patterns.

2

seq_along(vars) gives 1:length(vars), and is safer because it still behaves correctly if vars happens to be empty.

    age d.motor glucose 
  33.21    3.99    8.91 

The second pattern grows a data frame row by row with rbind(). This is slower, but it is easy to read, and Practical 3 uses it to collect the predictions of four different models.

1results <- NULL
for (v in vars) {
2  fit_v <- glm(tbi$d.unfav ~ tbi[[v]], family = binomial)
3  results <- rbind(results,
                   data.frame(variable = v,
                              coefficient = round(coef(fit_v)[[2]], 4),
                              aic = round(AIC(fit_v), 1)))
}
results

1

Start from an empty object.

2

Fit a univariable model for each predictor in turn.

3

Stack the new one-row data frame underneath what we have so far.

  variable coefficient    aic
1      age      0.0316 2806.0
2  d.motor     -0.6887 2622.7
3  glucose      0.1399 2704.6

5.2.3 The apply family

sapply() does the same as a for-loop, but in a single line, and returns the results directly.

1sapply(vars, \(v) mean(tbi[[v]], na.rm = TRUE)) |> round(2)

2sapply(tbi[, vars], class)

3lapply(vars, \(v) summary(tbi[[v]]))[[1]]

1

Apply the anonymous function to each element of vars and simplify the result into a named vector. Compare this with the six-line loop above.

2

Because a data frame is a list of columns, sapply() over a data frame applies the function to each column.

3

lapply() does the same, but always returns a list without simplifying it. This is useful when the results do not all have the same shape.

    age d.motor glucose 
  33.21    3.99    8.91 
      age   d.motor   glucose 
"integer" "integer" "numeric" 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  14.00   22.00   30.00   33.21   43.00   79.00 

Loop or sapply()?

Both do the same job. sapply() is convenient when each iteration produces one value and the iterations are independent of each other. A for-loop is easier to follow when each step depends on the previous one, or when the body of the loop is long.

5.2.4 Putting it together: a bootstrap

The code below is the pattern behind bootstrap validation, and a shortened version of what Practical 3 does to estimate a uniform shrinkage factor. It combines most of what we have covered so far.

1set.seed(1)

2dev <- tbi[!is.na(tbi$d.pupil), ]
fit_dev <- glm(d.unfav ~ age + factor(d.motor),
               data = dev, family = binomial)

3b <- 100
4slopes <- sapply(1:b, \(i) {

5  smp <- dev[sample(nrow(dev), replace = TRUE), ]

6  fit_b <- glm(d.unfav ~ age + factor(d.motor),
               data = smp, family = binomial)

7  lp <- predict(fit_b, newdata = dev, type = "link")

8  coef(glm(dev$d.unfav ~ lp, family = binomial))[[2]]
})

cat("Mean calibration slope over", b, "bootstraps:",
9    round(mean(slopes), 3), "\n")

1

Fix the seed so the whole procedure is reproducible.

2

Use only complete cases, so that every bootstrap sample can be fitted.

3

The number of bootstrap repetitions. Practical 3 uses 500.

4

For each of b repetitions, run the anonymous function and collect its single returned value.

5

Draw a sample of row indices with replacement, of the same size as the original data: one bootstrap sample.

6

Redevelop the model from scratch in that sample.

7

Apply the bootstrap model to the original data and obtain linear predictors.

8

Regress the original outcome on those linear predictors; the coefficient is the calibration slope.

9

Averaged over all repetitions, this is the shrinkage factor. A value below 1 indicates optimism, meaning that the model fits its own development sample better than it fits new data.

Mean calibration slope over 100 bootstraps: 0.989 

Exercise 7

1. Rewrite the means for-loop from above using sapply() in a single line.
2. Use sapply() to count the missing values in every column of tbi, and show only the columns that have any.
3. Adapt the bootstrap above to also record the intercept of each bootstrap calibration model. (Hint: return a vector of two values from the anonymous function and inspect the shape of the result.)
4. Increase b from 100 to 500 and see how much the mean slope moves. What does that tell you about how many bootstraps are enough?

Solution 7

# a.
sapply(vars, \(v) mean(tbi[[v]], na.rm = TRUE)) |> round(2)

    age d.motor glucose 
  33.21    3.99    8.91 

# b.
miss <- sapply(tbi, \(x) sum(is.na(x)))
miss[miss > 0]

 d.pupil  hypoxia hypotens  ctclass     tsah      edh cisterns    shift 
     123      249       59       25       96       38      257      252 
 glucose glucoset       ph   sodium  sodiumt       hb      hbt 
      64       64      457       62       62       24       24 

# c. returning two values per iteration gives a 2-by-b matrix
set.seed(1)
both <- sapply(1:100, \(i) {
  smp   <- dev[sample(nrow(dev), replace = TRUE), ]
  fit_b <- glm(d.unfav ~ age + factor(d.motor),
               data = smp, family = binomial)
  lp    <- predict(fit_b, newdata = dev, type = "link")
  cf    <- coef(glm(dev$d.unfav ~ lp, family = binomial))
  c(intercept = cf[[1]], slope = cf[[2]])
})

dim(both)

[1]   2 100

rowMeans(both) |> round(3)

intercept     slope 
   -0.007     0.989 

# d.
set.seed(1)
slopes_500 <- sapply(1:500, \(i) {
  smp   <- dev[sample(nrow(dev), replace = TRUE), ]
  fit_b <- glm(d.unfav ~ age + factor(d.motor),
               data = smp, family = binomial)
  lp    <- predict(fit_b, newdata = dev, type = "link")
  coef(glm(dev$d.unfav ~ lp, family = binomial))[[2]]
})
round(c(b100 = mean(slopes), b500 = mean(slopes_500)), 4)

  b100   b500 
0.9891 0.9847 

The estimate hardly moves, because we are taking an average. The mean stabilises well before the individual bootstrap estimates do, so a larger number of bootstraps mainly adds precision in the last decimal.

6 Part V: Plotting with {ggplot2}

Practicals 2 and 3 present their results graphically, using calibration plots, ROC curves, and effect plots, all of which are made with {ggplot2}. Such a plot always consists of the data, an aesthetic mapping from variables to visual properties, and one or more geometries that are drawn. Layers are added with +.

library(ggplot2)

1ggplot(data = tbi, aes(x = age)) +
2  geom_histogram(binwidth = 5,
                 fill = "#648FFF",
                 colour = "white") +
3  labs(x = "Age (yrs)", y = "Number of patients") +
4  theme_bw()

1

The data, and the mapping: age goes on the x-axis.

2

The geometry: draw it as a histogram.

3

Axis labels. labs() also takes title and subtitle.

4

A theme controls the overall appearance.

Adding a second variable is a matter of extending the mapping. The plot below is a first calibration plot, of the kind that Practical 2 builds:

1dev$preds <- fit_dev[["fitted.values"]]

2ggplot(dev, aes(x = preds, y = d.unfav)) +
3  geom_abline(linewidth = 1, colour = "black", alpha = 0.33) +
4  geom_point(alpha = 0.2) +
5  geom_smooth(method = "loess", formula = "y ~ x",
              colour = "#785EF0", fill = "#785EF0") +
6  coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
  labs(x = "Predicted probability", y = "Observed probability") +
  theme_bw()

1

Store the predicted probabilities as a column, so that ggplot() can find them.

2

Predicted probability on the x-axis, observed outcome on the y-axis.

3

The diagonal line of perfect calibration, drawn 67% transparent (alpha = 0.33).

4

One point per patient. The transparency keeps overlapping points readable.

5

A smoothed curve through the points: this is the calibration curve.

6

Restrict both axes to the range of a probability.

Two further features that we need are mapping a variable to colour and splitting a plot into panels.

1ggplot(dev, aes(x = age, y = preds, colour = factor(d.unfav))) +
  geom_point(alpha = 0.4) +
2  scale_colour_manual(
    values = c("#648FFF", "#FFB000"),
    labels = c("Favourable", "Unfavourable")
  ) +
3  facet_wrap(~ trial) +
  labs(x = "Age (yrs)", y = "Predicted probability") +
  theme_bw() +
4  theme(legend.position = "bottom", legend.title = element_blank())

1

Mapping d.unfav to colour splits the points into two coloured groups. The variable has to be a factor, otherwise {ggplot2} draws a continuous colour scale.

2

Set the colours and the legend labels manually.

3

facet_wrap() draws one panel per level of a variable, using a one-sided formula.

4

Some fine-tuning, by moving the legend and removing its title.

Exercise 8

1. Draw a boxplot of age by d.pupil. (Hint: geom_boxplot().)
2. Draw a histogram of the predicted probabilities in dev, coloured by the observed outcome.
3. Plot predicted probability against age, with a LOESS smoother, faceted by cause.

Solution 8

# a.
ggplot(dev, aes(x = d.pupil, y = age)) +
  geom_boxplot(fill = "#648FFF", alpha = 0.5) +
  labs(x = "Pupillary reactivity", y = "Age (yrs)") +
  theme_bw()

# b.
ggplot(dev, aes(x = preds, fill = factor(d.unfav))) +
  geom_histogram(binwidth = 0.02, position = "identity", alpha = 0.6) +
  scale_fill_manual(values = c("#648FFF", "#FFB000"),
                    labels = c("Favourable", "Unfavourable")) +
  labs(x = "Predicted probability", y = "Count") +
  theme_bw() +
  theme(legend.position = "bottom", legend.title = element_blank())

# c.
ggplot(dev, aes(x = age, y = preds)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess", formula = "y ~ x", colour = "#785EF0") +
  facet_wrap(~ cause) +
  labs(x = "Age (yrs)", y = "Predicted probability") +
  theme_bw()

7 Part VI: Reporting results

Model output usually needs some tidying before it goes into a report. The functions below are used throughout Practicals 2 and 3.

auc_value <- 0.7654321

1round(auc_value, 3)
2format(round(0.7, 3), nsmall = 3)
3paste0("AUC = ", round(auc_value, 3))
4cat("AUC:", round(auc_value, 3), "\n")
5cat("Intercept:", 0.001, "\nSlope:", 1.002, "\n")
6sprintf("AUC = %.3f (95%% CI %.3f-%.3f)", 0.765, 0.74, 0.79)

1

Round to three decimals.

2

round() drops trailing zeros; format(..., nsmall = 3) forces them back, so that columns line up.

3

paste0() pastes strings together without a separator, whereas paste() uses a space.

4

cat() prints without quotes or index markers, which is what we want for a readable message. The \n stands for a line break.

5

Several values in one call, with a line break in the middle.

6

sprintf() gives precise control over formatting. %.3f means “a number with three decimals”; %% prints a literal percent sign.

[1] 0.765
[1] "0.700"
[1] "AUC = 0.765"
AUC: 0.765 
Intercept: 0.001 
Slope: 1.002 
[1] "AUC = 0.765 (95% CI 0.740-0.790)"

print() shows the object as R stores it, including quotes and [1] markers, whereas cat() prints the text as it should be read. In general, we use cat() for messages and print() for objects.

8 Capstone: from raw data to a fitted model

This last section combines the previous parts into one short analysis, and stops where Practical 2 begins. Work through it line by line and check that every step is clear.

# ---- 1. Start clean and set a seed -------------------------------
1rm(list = ls())
set.seed(1)

tbi <- rio::import("I:/onderwijs/PREDICT (voormalig PINT)/2026 Najaar/R practica/TBI_data.txt")

# ---- 2. Select outcome and predictors ----------------------------
2outcome    <- "d.unfav"
predictors <- c("age", "d.motor", "d.pupil", "cause")

3analysis <- tbi[, c(outcome, predictors)]
str(analysis)

1

Clearing the workspace prevents an object from an earlier session from affecting the results.

2

Storing the names in variables makes the script easier to adjust later. This works because [[ and [ accept character vectors.

3

Keep only what the analysis needs.

'data.frame':   2159 obs. of  5 variables:
 $ d.unfav: int  0 0 0 0 0 0 0 0 0 0 ...
 $ age    : int  14 14 14 14 14 14 14 14 14 14 ...
 $ d.motor: int  5 4 4 4 5 3 5 5 4 5 ...
 $ d.pupil: chr  "both reactive" "both reactive" "both reactive" "both reactive" ...
 $ cause  : chr  "Motorbike" "Motorbike" "domestic/fall" "Motorbike" ...

# ---- 3. Inspect the outcome --------------------------------------
1table(analysis[[outcome]])
round(proportions(table(analysis[[outcome]])) * 100, 1)

# ---- 4. Inspect missingness --------------------------------------
2colSums(is.na(analysis))

1

The number of events, which largely determines the risk of overfitting.

2

Only d.pupil is incomplete. Practical 2 imputes it; here we simply note it.

   0    1 
1308  851 

   0    1 
60.6 39.4 
d.unfav     age d.motor d.pupil   cause 
      0       0       0     123       0 

# ---- 5. Declare categorical variables as factors ------------------
1analysis[["d.motor"]] <- as.factor(analysis[["d.motor"]])
analysis[["d.pupil"]] <- factor(
  analysis[["d.pupil"]],
  levels = c("no reactive pupils", "one reactive", "both reactive"),
2  labels = 0:2)
analysis[["cause"]] <- factor(analysis[["cause"]])

3sapply(analysis, class)

1

Numeric codes must be declared as factors, or they will be modelled as continuous.

2

Choosing the level order explicitly fixes both the reference category and the meaning of the labels.

3

Check that every variable now has the type that we intended.

  d.unfav       age   d.motor   d.pupil     cause 
"integer" "integer"  "factor"  "factor"  "factor" 

# ---- 6. Fit the model --------------------------------------------
1f <- as.formula(paste(outcome, "~", paste(predictors, collapse = " + ")))
f

2fit <- glm(f, data = analysis, family = binomial)
summary(fit)

1

Build the formula from the variable names we stored in step 2. paste(..., collapse = " + ") joins the predictors with plus signs, and as.formula() turns the resulting string into a formula object.

2

Fit the logistic regression model.

d.unfav ~ age + d.motor + d.pupil + cause

Call:
glm(formula = f, family = binomial, data = analysis)

Coefficients:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)                -0.347931   0.664542  -0.524   0.6006    
age                         0.037503   0.004037   9.289  < 2e-16 ***
d.motor2                    0.648468   0.621264   1.044   0.2966    
d.motor3                    0.057519   0.614843   0.094   0.9255    
d.motor4                   -0.660139   0.609902  -1.082   0.2791    
d.motor5                   -1.304401   0.610827  -2.135   0.0327 *  
d.motor6                   -1.628588   0.680168  -2.394   0.0166 *  
d.pupil1                   -0.733449   0.185156  -3.961 7.46e-05 ***
d.pupil2                   -1.285761   0.148476  -8.660  < 2e-16 ***
causedomestic/fall          0.326742   0.247825   1.318   0.1874    
causeMotorbike              0.162727   0.245626   0.662   0.5077    
causeother                  0.435268   0.245729   1.771   0.0765 .  
causeRoad traffic accident  0.183664   0.231363   0.794   0.4273    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2726.7  on 2035  degrees of freedom
Residual deviance: 2251.5  on 2023  degrees of freedom
  (123 observations deleted due to missingness)
AIC: 2277.5

Number of Fisher Scoring iterations: 4

# ---- 7. Reach into the fit ---------------------------------------
1length(coef(fit))
2length(fit[["y"]])

3round(exp(coef(fit)[["age"]]), 3)

analysis[["preds"]] <- NA
4analysis[["preds"]][!is.na(tbi$d.pupil)] <- fit[["fitted.values"]]

cat("Number of coefficients:", length(coef(fit)),
    "\nPatients used:", length(fit[["y"]]),
    "\nMean predicted risk:", round(mean(fit[["fitted.values"]]), 3),
5    "\nObserved risk:", round(mean(fit[["y"]]), 3), "\n")

1

Thirteen coefficients, because three of the four predictors are categorical. This is the number Practical 2 asks you to verify.

2

Fewer than 2,159 patients: the ones with a missing d.pupil were silently dropped.

3

The odds ratio for one extra year of age, obtained by exponentiating the log odds ratio.

4

Write the predictions back into the data frame, taking care to place them in the rows the model actually used.

5

A short and readable summary. The mean predicted risk and the observed risk agree, as they should for a correctly fitted logistic model on its own development data.

[1] 13
[1] 2036
[1] 1.038
Number of coefficients: 13 
Patients used: 2036 
Mean predicted risk: 0.392 
Observed risk: 0.392 

Where Practical 2 picks up

We now have a fitted model, together with its coefficients, predicted probabilities, and linear predictors, and we know how to get each of these out of the fit object. Practical 2 starts at this point and continues with the next questions. Is age linearly related to the outcome? What should we do about the patients that were dropped? And how well does this model discriminate and calibrate?

9 Getting unstuck

Errors are a normal part of working in R. The table below lists a few common ones, together with what to check first.

A first-aid checklist for common R errors

Problem	First thing to try
“could not find function”	The package is not loaded. Use library() or pkg::fun().
“object not found”	A typo, or the object was never created. Check ls().
“undefined columns selected”	You are subsetting a column that does not exist. Check names().
Unexpected NAs	A coercion failed, or na.rm = TRUE is missing.
The result has the wrong shape	Check with class(), dim(), and str(). Usually [ was used where [[ was needed.
A factor behaves oddly	Check levels(). Convert with as.numeric(as.character(x)).
No idea what an object is	str(x, max.level = 1).

Apart from that, ?glm opens the help page of a function and ??calibration searches all help pages. The Examples section at the bottom of a help page is often the quickest way to understand an unfamiliar function.

10 Final words

We have now covered the building blocks of R that we need for prediction modelling: packages and reproducibility, vectors and types, subsetting with [, [[, and $, factors and reference categories, lists and the structure of a model fit object, formulas, data wrangling, functions with default arguments and conditionals, loops and the apply family, a first bootstrap, plotting with {ggplot2}, and formatting results.

None of this was statistics, but it is what the next two practicals build on. From Practical 2 onwards, R itself should no longer be the difficult part.

Before the next practical, make sure that {pacman} is installed, and have a look at the tutorial on how to create a prediction model: https://github.com/CHMMaas/TutorialPredictionModel

11 Postscript

This practical was developed on:

R version 4.4.0 (2024-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)

Matrix products: default


locale:
[1] LC_COLLATE=Dutch_Netherlands.utf8  LC_CTYPE=Dutch_Netherlands.utf8   
[3] LC_MONETARY=Dutch_Netherlands.utf8 LC_NUMERIC=C                      
[5] LC_TIME=Dutch_Netherlands.utf8    

time zone: Europe/Amsterdam
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggplot2_4.0.3 dplyr_1.1.4   stringr_1.6.0

loaded via a namespace (and not attached):
 [1] Matrix_1.7-0       gtable_0.3.6       jsonlite_2.0.0     compiler_4.4.0    
 [5] tidyselect_1.2.1   splines_4.4.0      scales_1.4.0       yaml_2.3.12       
 [9] fastmap_1.2.0      lattice_0.22-6     mime_0.13          R6_2.6.1          
[13] labeling_0.4.3     generics_0.1.4     knitr_1.51         htmlwidgets_1.6.4 
[17] tibble_3.2.1       pillar_1.11.1      RColorBrewer_1.1-3 R.utils_2.13.0    
[21] rlang_1.2.0        stringi_1.8.7      xfun_0.57          S7_0.2.2          
[25] otel_0.2.0         cli_3.6.5          mgcv_1.9-1         withr_3.0.2       
[29] magrittr_2.0.3     rio_1.2.4          digest_0.6.37      grid_4.4.0        
[33] rstudioapi_0.18.0  nlme_3.1-164       lifecycle_1.0.5    R.methodsS3_1.8.2 
[37] R.oo_1.27.1        vctrs_0.7.3        evaluate_1.0.5     glue_1.8.0        
[41] data.table_1.17.0  farver_2.1.2       rmarkdown_2.31     tools_4.4.0       
[45] pkgconfig_2.0.3    htmltools_0.5.8.1