1.6.1 Missing values

Real data always has gaps. A student skips a question, or a value is lost. R marks a missing value with NA, short for not available. It is neither zero nor empty text; it means “we do not know”:

sleep_with_gap <- c(6.5, NA, 5.5)
mean(sleep_with_gap)

[1] NA

If one value is unknown, the average is unknown too, so R returns NA rather than guessing. You will see shortly how to tell R to leave missing values out.

1.7.1 Getting help

Every function has a help page. Type a question mark before its name in the Console:

?mean

The page opens in the Help pane. Help pages are written for experienced users, so they can look dense at first. Start with three parts: Usage (how to call the function), Arguments (what each input means), and Examples at the bottom, which you can copy and run.

When a search engine or an AI assistant suggests code, treat it like advice from a knowledgeable stranger: often right, sometimes wrong, and always worth checking. Chapter 18 shows how to use AI tools well.

1.12.1 Summary

• Every analysis is a chain of decisions. Writing it as code records each one, so the analysis can be checked and repeated: it becomes reproducible.
• R is a free language made for working with data. You write it in RStudio: the Source pane holds your scripts, the Console runs the code, the Environment pane lists your objects, and the fourth pane shows files, plots, packages, and help.
• <- stores a value in an object. c() combines several values.
• Every value has a type. The most common are numeric, character (in quotes), and logical (TRUE or FALSE). NA marks a missing value.
• Functions take arguments inside brackets. ? opens a function’s help page.
• Write your code in scripts, keep each piece of work in an RStudio Project, install packages once with install.packages(), and load them each session with library().

1.12.2 Key terms

R, RStudio, Console, script, comment, object, assignment, vector, data type, missing value (NA), function, argument, package, data frame, working directory, RStudio Project, reproducible research.

2.2.1 Operations on a whole vector

Most things you do to a vector happen to every value at once. Converting the sleep values from hours into minutes needs only one multiplication:

sleep * 60

[1] 390 420 330 480 360

Comparisons work the same way, and give one TRUE or FALSE for each value. The nights shorter than the recommended 7 hours are marked by:

sleep < 7

[1]  TRUE FALSE  TRUE FALSE  TRUE

Because R counts TRUE as 1 and FALSE as 0, sum() counts the short nights, and mean() gives their share:

sum(sleep < 7)

[1] 3

mean(sleep < 7)

[1] 0.6

Three of the five nights, or 60%, were short. The same line of code would work just as well on 2,000 nights.

2.2.2 Selecting values

Square brackets pick values from a vector by their position. Positions start at 1:

sleep[1]         # the first night

[1] 6.5

sleep[c(1, 3)]   # the first and third nights

[1] 6.5 5.5

sleep[-2]        # every night except the second

[1] 6.5 5.5 8.0 6.0

Values can also be picked with a condition, by putting a TRUE/FALSE vector inside the brackets. R keeps the values where the condition is TRUE:

sleep[sleep < 7]

[1] 6.5 5.5 6.0

The line reads aloud as “sleep, where sleep is less than 7”. This way of selecting, called logical indexing, is one of the most useful ideas in R, and you will use it constantly.

2.2.3 Changing values

Brackets on the left of the arrow change values. Suppose the second night turns out to have been 7.5 hours, not 7:

sleep[2] <- 7.5
sleep

[1] 6.5 7.5 5.5 8.0 6.0

2.3.1 Kinds of measurement and R’s types

The choice between a factor and a number is not a matter of convenience; it follows from the kind of measurement that produced the values. Categories without an order, such as faculty, are nominal. Categories with an order but uneven or unknown distances between them, such as none, part-time, and full-time employment, are ordinal. Measurements on a scale with equal distances, such as hours of sleep or a wellbeing index, are numeric. Table 2.1 shows how each kind is stored in R. Chapter 5 explains these levels of measurement in full, and why they decide which summaries are meaningful.

2.4.1 Selecting columns and rows

A single column is picked with $, as in Chapter 1:

pilot$sleep

[1] 6.5 7.5 5.5 8.0

Square brackets work on data frames too, with two positions separated by a comma: rows first, then columns. An empty position means “all”:

pilot[2, ]                        # the second row, all columns

  student    faculty sleep invited
2      S2 Humanities   7.5   FALSE

pilot[, c("student", "sleep")]    # all rows, two columns

  student sleep
1      S1   6.5
2      S2   7.5
3      S3   5.5
4      S4   8.0

pilot[pilot$sleep < 7, ]          # rows where sleep is under 7

  student   faculty sleep invited
1      S1 Education   6.5    TRUE
3      S3 Education   5.5    TRUE

The last line combines a data frame with logical indexing: “pilot, the rows where sleep is under 7, all columns”. It is the R version of SPSS’s Select Cases.

2.4.2 Adding columns

Assigning to a new column name adds a column:

pilot$short_sleep <- pilot$sleep < 7
pilot

  student         faculty sleep invited short_sleep
1      S1       Education   6.5    TRUE        TRUE
2      S2      Humanities   7.5   FALSE       FALSE
3      S3       Education   5.5    TRUE        TRUE
4      S4 Health Sciences   8.0   FALSE       FALSE

2.4.3 The student wellbeing data

The students data frame from the data2thesis package has the same structure as pilot, only larger:

str(students)

'data.frame':   600 obs. of  14 variables:
 $ student_id         : chr  "S0001" "S0002" "S0003" "S0004" ...
 $ supervisor_id      : chr  "SUP088" "SUP116" "SUP040" "SUP112" ...
 $ age                : int  31 28 33 34 35 36 NA 25 29 26 ...
 $ gender             : chr  "Female" "Female" "Male" "Female" ...
 $ faculty            : chr  "Health Sciences" "Education" "Education" "Health Sciences" ...
 $ programme          : chr  "Master's" "PhD" "PhD" "Master's" ...
 $ study_mode         : chr  "Full-time" "Full-time" "Full-time" "Part-time" ...
 $ employment         : chr  "None" "Part-time job" "None" "Full-time job" ...
 $ has_children       : chr  "No" "No" "No" "No" ...
 $ lives_away         : chr  "Yes" "Yes" "No" "No" ...
 $ financial_worry    : int  2 2 4 2 2 5 3 NA 2 2 ...
 $ workshop           : chr  "Not invited" "Invited" "Invited" "Invited" ...
 $ workshop_sessions  : int  0 0 5 6 2 0 4 0 0 0 ...
 $ considering_dropout: chr  "No" "No" "No" "No" ...

In this listing, chr means character, int whole numbers, and num numbers with decimals. The categories, such as faculty, arrive as text, and they are converted into factors when they are needed as categories:

students$faculty <- factor(students$faculty)
levels(students$faculty)

[1] "Education"        "Health Sciences"  "Humanities"       "Natural Sciences"
[5] "Social Sciences" 

Selection works as before. The number of PhD students in Health Sciences, for example, is found with two conditions:

hs_phd <- students[students$faculty == "Health Sciences" & students$programme == "PhD", ]
nrow(hs_phd)

[1] 44

Two symbols are new here. The double equals sign, ==, means “is equal to”; a single = is used for arguments, so comparisons need two. The ampersand, &, means “and”: both conditions must be TRUE. Its partner | means “or”.

2.5.1 Matrices

A matrix is a grid of values of one type, usually numbers. Matrices appear most often as results. The function cor(), for example, calculates the correlations between several variables and returns them as a matrix. The correlations below use the first semester only, so that each student counts once:

first_semester <- semesters[semesters$semester == 1, ]
vars <- first_semester[, c("gpa", "sleep_hours", "study_hours", "wellbeing")]
correlations <- cor(vars, use = "complete.obs")
round(correlations, 2)

             gpa sleep_hours study_hours wellbeing
gpa         1.00        0.25        0.07      0.31
sleep_hours 0.25        1.00       -0.56      0.46
study_hours 0.07       -0.56        1.00     -0.24
wellbeing   0.31        0.46       -0.24      1.00

Every variable correlates perfectly with itself, hence the 1s on the diagonal. Chapter 6 explains how to read correlations. For now, notice that values are picked from a matrix exactly as from a data frame, with [row, column]:

correlations["sleep_hours", "wellbeing"]

[1] 0.4602433

2.5.2 Lists

A list is a container that can hold anything: vectors of different lengths, data frames, even other lists. Its parts are usually named:

student <- list(
  id        = "S0001",
  programme = "Master's",
  sleep     = c(7.2, 8.4, 7.8, 7.2)
)
student$sleep

[1] 7.2 8.4 7.8 7.2

Lists are rarely built by hand, but R’s statistical functions return their results as lists. A preview of a test from Chapter 7, which examines whether students sleep 7 hours on average in their first semester, shows this:

result <- t.test(first_semester$sleep_hours, mu = 7)
names(result)

 [1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"   
 [6] "null.value"  "stderr"      "alternative" "method"      "data.name"  

result$estimate

mean of x 
 6.483502 

The function names() lists everything the test calculated, and $ picks out one part, here the students’ average sleep. This is how a number from a statistical test is retrieved for a report.

2.6.1 CSV files

A CSV file (comma-separated values) is a plain text table, where commas separate the columns. It is the most common format for sharing data, and R reads it with read.csv():

path <- data2thesis_example("semesters.csv")
semester_file <- read.csv(path)
head(semester_file, 3)

  student_id semester  gpa sleep_hours study_hours exercise_days caffeine_mg
1      S0001        1 3.02         7.2          26             5         130
2      S0001        2 3.27         8.4          17             3          80
3      S0001        3 3.17         7.8          25             7         105
  supervisor_meetings wellbeing
1                   0        74
2                   1        69
3                   3        71

In your own project, this would simply be read.csv("semesters.csv"), or read.csv(here::here("data", "semesters.csv")) with a data folder.

2.6.2 Excel files

The readxl package reads Excel files. Install it once with install.packages("readxl"):

library(readxl)
raw <- read_excel(data2thesis_example("wellbeing_raw.xlsx"))
dim(raw)

[1] 615  64

This is the export from the study’s survey tool, and it shows what real survey exports look like. The function dim() gives its size: 615 rows (more than the 600 students) and 64 columns. The first rows are:

head(raw[, 1:6])

# A tibble: 6 × 6
  `Response ID` `Q1_Student ID` Q2_Age Q3_Gender Q4_Faculty      Q5_Programme
  <chr>         <chr>           <chr>  <chr>     <chr>           <chr>       
1 R0001         TEST            99     Female    <NA>            <NA>        
2 R0002         test            <NA>   <NA>      <NA>            <NA>        
3 R0003         TEST2           30     Male      <NA>            <NA>        
4 R0004         S0191           27     male      Humanities      Masters     
5 R0005         S0015           30     Male      Health Sciences Master's    
6 R0006         S0474           38     F         health sciences Masters     

The printout looks a little different from before, because read_excel() returns a tibble: a modern kind of data frame that prints more compactly and shows each column’s type under its name. It can be used exactly like a data frame.

Several problems are already visible. There are test responses (TEST), odd codes such as 99 for missing age, and the same answer spelled several ways (Male, male). Every column shown is <chr>, text, even the ages, because a few entries are not plain numbers. Chapter 3 is devoted to cleaning this file; for now, it is enough to be able to open it.

If a workbook has several sheets, choose one with the sheet argument, for example read_excel("file.xlsx", sheet = "Semester 2").

2.6.3 SPSS files

The haven package reads SPSS files (.sav), and also Stata and SAS files:

library(haven)
spss <- read_sav(data2thesis_example("wellbeing.sav"))
spss[1:3, c("student_id", "gender", "faculty", "stress_1")]

# A tibble: 3 × 4
  student_id gender     faculty             stress_1   
  <chr>      <dbl+lbl>  <dbl+lbl>           <dbl+lbl>  
1 S0001      1 [Female] 2 [Health Sciences] 3 [Neutral]
2 S0002      1 [Female] 1 [Education]       4 [Agree]  
3 S0003      2 [Male]   1 [Education]       4 [Agree]  

SPSS files carry labels, and haven keeps them. Each value is stored as a number, with its label shown in square brackets: 2 [Health Sciences]. The question behind each variable is kept too:

attr(spss$stress_1, "label")

[1] "I feel unable to control important things in my studies."

To analyse a labelled variable as categories, turn it into a factor with as_factor(), which uses the labels as levels:

table(as_factor(spss$faculty))

       Education  Health Sciences       Humanities Natural Sciences 
             148              154               95               87 
 Social Sciences 
             116 

This is one of the great conveniences for anyone moving from SPSS: the work put into labelling the data is not lost.

2.12.1 Summary

• Data is the record of measurements: a table of cases (rows) and variables (columns). The unit of analysis is what one row represents, and the same study can have more than one.
• A vector holds the values of one variable, all of one type. Mixing types turns everything into text.
• Most operations work on every value of a vector at once. [ ] picks values by position or by a condition (logical indexing).
• A factor holds categories with a fixed set of levels, whose order you can set. The kind of measurement decides whether a variable should be a factor or a number.
• A data frame is a table: one row per case, one column (a vector) per variable. Select with $ or [rows, columns], and use str() to see its structure.
• A matrix is a grid of one type, often a result such as a correlation matrix. A list can hold anything; statistical results are lists.
• read.csv(), read_excel() (readxl), and read_sav() (haven) import data. haven keeps SPSS labels, and as_factor() turns them into factors. write.csv(), write_xlsx(), and write_sav() export data.
• A codebook records what every variable means and how it was measured; SPSS labels provide much of it.
• Error messages say what R could not do. Read them, and check spelling and capitals first.

2.12.2 Key terms

Case, variable, data matrix, unit of analysis, data structure, vector, type conversion, logical indexing, factor, level, data frame, tibble, matrix, list, CSV file, variable label, value label, codebook, error message.

3.3.1 Chaining steps with the pipe

Analysis is a series of steps: take the data, keep some rows, calculate something, round the result. In Chapter 1, such steps were written by putting functions inside each other:

sleep <- c(6.5, 7, 5.5, 8, 6)
round(mean(sleep), 1)

[1] 6.6

This reads from the inside out, which gets hard to follow as the steps multiply. The pipe, written |>, writes the steps in the order they happen. It takes the result on its left and passes it as the first argument to the function on its right:

sleep |> mean() |> round(1)

[1] 6.6

The pipe reads as “and then”: take sleep, and then take the mean, and then round it to one decimal place. In RStudio, Ctrl+Shift+M (Cmd+Shift+M on a Mac) types the pipe for you.

3.4.1 Choosing the cases you need

Most analyses use only some of the cases: the students of one faculty, the first semester, the participants who completed the study. The verb filter() keeps the rows that meet a condition:

pilot |> filter(sleep < 7)

# A tibble: 3 × 4
  student faculty    sleep stress
  <chr>   <chr>      <dbl>  <dbl>
1 S1      Education    6.5    3.2
2 S3      Education    5.5    4.5
3 S5      Humanities   6      3.9

Several conditions separated by commas must all be true:

pilot |> filter(sleep < 7, faculty == "Education")

# A tibble: 2 × 4
  student faculty   sleep stress
  <chr>   <chr>     <dbl>  <dbl>
1 S1      Education   6.5    3.2
2 S3      Education   5.5    4.5

To keep rows matching any of several values, use %in%:

pilot |> filter(faculty %in% c("Education", "Humanities"))

# A tibble: 4 × 4
  student faculty    sleep stress
  <chr>   <chr>      <dbl>  <dbl>
1 S1      Education    6.5    3.2
2 S2      Humanities   7.5    2.1
3 S3      Education    5.5    4.5
4 S5      Humanities   6      3.9

Every filter is a decision about the sample, and should be reported: an analysis of “students who completed all four semesters” answers a different question from an analysis of all students.

3.4.2 Choosing variables

Large datasets have far more variables than any one analysis needs. The verb select() keeps the columns you name, in the order you name them, and a minus sign drops a column instead:

pilot |> select(student, sleep)

# A tibble: 5 × 2
  student sleep
  <chr>   <dbl>
1 S1        6.5
2 S2        7.5
3 S3        5.5
4 S4        8  
5 S5        6  

pilot |> select(-faculty)

# A tibble: 5 × 3
  student sleep stress
  <chr>   <dbl>  <dbl>
1 S1        6.5    3.2
2 S2        7.5    2.1
3 S3        5.5    4.5
4 S4        8      1.8
5 S5        6      3.9

3.4.3 Sorting cases

Sorting helps when looking at data, for example to see the most extreme values first. The verb arrange() sorts rows by one or more columns, and desc() sorts from largest to smallest:

pilot |> arrange(desc(stress))

# A tibble: 5 × 4
  student faculty         sleep stress
  <chr>   <chr>           <dbl>  <dbl>
1 S3      Education         5.5    4.5
2 S5      Humanities        6      3.9
3 S1      Education         6.5    3.2
4 S2      Humanities        7.5    2.1
5 S4      Health Sciences   8      1.8

3.4.4 Creating new variables

Many variables used in analysis are calculated from others: a total score, a conversion to other units, a category derived from a number. The verb mutate() adds new columns calculated from existing ones, or changes existing ones:

pilot |> mutate(
  sleep_minutes = sleep * 60,
  short_sleep   = sleep < 7
)

# A tibble: 5 × 6
  student faculty         sleep stress sleep_minutes short_sleep
  <chr>   <chr>           <dbl>  <dbl>         <dbl> <lgl>      
1 S1      Education         6.5    3.2           390 TRUE       
2 S2      Humanities        7.5    2.1           450 FALSE      
3 S3      Education         5.5    4.5           330 TRUE       
4 S4      Health Sciences   8      1.8           480 FALSE      
5 S5      Humanities        6      3.9           360 TRUE       

To sort values into categories, case_when() checks conditions in order and uses the first one that is true, and .default catches everything else:

pilot |> mutate(
  stress_level = case_when(
    stress >= 4 ~ "High",
    stress >= 2.5 ~ "Medium",
    .default = "Low"
  )
)

# A tibble: 5 × 5
  student faculty         sleep stress stress_level
  <chr>   <chr>           <dbl>  <dbl> <chr>       
1 S1      Education         6.5    3.2 Medium      
2 S2      Humanities        7.5    2.1 Low         
3 S3      Education         5.5    4.5 High        
4 S4      Health Sciences   8      1.8 Low         
5 S5      Humanities        6      3.9 Medium      

Turning a number into categories in this way loses information, since a stress score of 3.9 and one of 2.6 become the same “Medium”. It is useful for description, but analyses should normally keep the original number.

3.4.5 Summarising, overall and by group

Summaries reduce many values to a few numbers. The verb summarise() reduces a table to one row of summaries, and n() counts the rows:

pilot |> summarise(
  mean_sleep = mean(sleep),
  students   = n()
)

# A tibble: 1 × 2
  mean_sleep students
       <dbl>    <int>
1        6.7        5

The real power comes with groups. The .by argument calculates the summary separately for each group:

pilot |> summarise(
  mean_sleep = mean(sleep),
  students   = n(),
  .by = faculty
)

# A tibble: 3 × 3
  faculty         mean_sleep students
  <chr>                <dbl>    <int>
1 Education             6           2
2 Humanities            6.75        2
3 Health Sciences       8           1

For simple counts, count() is a shortcut:

pilot |> count(faculty)

# A tibble: 3 × 2
  faculty             n
  <chr>           <int>
1 Education           2
2 Health Sciences     1
3 Humanities          2

3.4.6 Combining the steps

Each verb is simple, but chained together they answer real questions. One such question concerns how average wellbeing changed over the four semesters, separately for students who were and were not invited to the workshop. The semester records are in semesters, and the workshop group is in students, so the workshop group is first added to the semester records. (The left_join() in the first line is explained in Section 3.6.)

semesters |>
  left_join(students |> select(student_id, workshop), join_by(student_id)) |>
  summarise(
    mean_wellbeing = mean(wellbeing),
    students       = n(),
    .by = c(semester, workshop)
  ) |>
  arrange(semester, workshop)

  semester    workshop mean_wellbeing students
1        1     Invited       60.64667      300
2        1 Not invited       60.26667      300
3        2     Invited       65.47000      300
4        2 Not invited       60.22000      300
5        3     Invited       63.22261      283
6        3 Not invited       59.22857      280
7        4     Invited       60.93286      283
8        4 Not invited       58.42500      280

Read aloud, the code says: take the semester records, and then add each student’s workshop group, and then calculate the average wellbeing and the number of students for each semester and workshop group, and then sort the result. In semester 1, before the workshop, the two groups are almost level. From semester 2, the invited students are ahead. Chapter 7 tests whether that difference is larger than chance.

The number of students also drops in semesters 3 and 4. Some students left the study, which is worth remembering: it returns later in this chapter and in Chapter 6.

3.7.1 Inspecting the export

Cleaning starts with looking, because problems that have not been seen cannot be fixed:

library(readxl)
raw <- read_excel(data2thesis_example("wellbeing_raw.xlsx"))
dim(raw)

[1] 615  64

The file has 615 rows, although there are only 600 students. A few lines of code show what needs fixing:

raw |> count(Q3_Gender)

# A tibble: 7 × 2
  Q3_Gender     n
  <chr>     <int>
1 F            28
2 Female      233
3 M            46
4 Male        173
5 female       56
6 male         78
7 <NA>          1

raw |> filter(str_detect(str_to_upper(`Q1_Student ID`), "TEST")) |> select(1:4)

# A tibble: 3 × 4
  `Response ID` `Q1_Student ID` Q2_Age Q3_Gender
  <chr>         <chr>           <chr>  <chr>    
1 R0001         TEST            99     Female   
2 R0002         test            <NA>   <NA>     
3 R0003         TEST2           30     Male     

Gender is spelled six ways (F, female, Female, and so on), plus a blank from a test response, and there are test responses from when the survey was set up. (The raw file also contains "Female " with a trailing space; read_excel() removes such spaces automatically, but other import functions do not, which is why the cleaning below still uses str_trim().) Column names such as Q1_Student ID contain spaces, so they must be written inside backticks (`) in code. And as Chapter 2 showed, many numeric columns arrived as text.

3.7.2 Test responses and duplicates

Test responses are not data about students, and duplicated submissions would count some students twice, so both must go before anything else. Test responses have a student ID starting with “TEST”, in upper or lower case. The function str_to_upper() makes the case irrelevant, and str_detect() checks for the pattern. The function distinct() then removes rows that are exact copies of another row. The response ID is different for each submission, even for duplicates, so it is dropped first:

responses <- raw |>
  filter(!str_detect(str_to_upper(`Q1_Student ID`), "^TEST")) |>
  select(-`Response ID`) |>
  distinct()
nrow(responses)

[1] 600

In the pattern, ! means “not” and ^ means “at the start”. Exactly 600 rows are left: one per student.

3.7.3 Usable variable names

Names such as Q10_How worried are you about money? (1-5) are awkward to type and easy to get wrong. Short, consistent names make the code readable and reduce errors. The function rename() changes names, new name first, and the 22 questionnaire items, Q11_1 to Q11_22, are renamed in one go with rename_with(), using the scale each item belongs to:

item_names <- c(paste0("stress_", 1:6), paste0("burnout_", 1:6),
                paste0("support_", 1:6), paste0("satisfaction_", 1:4))

responses <- responses |>
  rename(
    student_id          = `Q1_Student ID`,
    age                 = Q2_Age,
    gender              = Q3_Gender,
    faculty             = Q4_Faculty,
    programme           = Q5_Programme,
    study_mode          = `Q6_Study mode`,
    employment          = `Q7_Paid work`,
    has_children        = Q8_Children,
    lives_away          = `Q9_Moved away from family`,
    financial_worry     = `Q10_How worried are you about money? (1-5)`,
    workshop            = `Workshop group`,
    workshop_sessions   = `Workshop sessions attended`,
    considering_dropout = `Y1_Considered leaving?`
  ) |>
  rename_with(~ item_names, .cols = Q11_1:Q11_22)

The call paste0("stress_", 1:6) builds the names stress_1 to stress_6, which saves typing them all. The questionnaire’s wording is not lost: it belongs in the codebook (Chapter 2).

3.7.4 Consistent categories

If “F”, “female”, and “Female” stay separate, every table and every test treats them as three different groups. Categories are fixed by looking for a pattern rather than listing every spelling, after str_to_lower() and str_trim() have removed differences in case and stray spaces:

responses <- responses |>
  mutate(
    gender = if_else(str_starts(str_to_lower(str_trim(gender)), "f"), "Female", "Male"),
    faculty = case_when(
      str_detect(str_to_lower(faculty), "educ")    ~ "Education",
      str_detect(str_to_lower(faculty), "health")  ~ "Health Sciences",
      str_detect(str_to_lower(faculty), "humanit") ~ "Humanities",
      str_detect(str_to_lower(faculty), "social")  ~ "Social Sciences",
      str_detect(str_to_lower(faculty), "natural|^sciences$") ~ "Natural Sciences"
    ),
    programme  = if_else(str_detect(str_to_lower(programme), "ph"), "PhD", "Master's"),
    study_mode = if_else(str_detect(str_to_lower(study_mode), "part"), "Part-time", "Full-time"),
    employment = if_else(str_to_lower(employment) %in% c("none", "no job"), "None", employment),
    across(c(has_children, lives_away, considering_dropout),
           ~ if_else(str_starts(str_to_lower(.x), "y"), "Yes", "No"))
  )
responses |> count(gender)

# A tibble: 2 × 2
  gender     n
  <chr>  <int>
1 Female   312
2 Male     288

responses |> count(faculty)

# A tibble: 5 × 2
  faculty              n
  <chr>            <int>
1 Education          148
2 Health Sciences    154
3 Humanities          95
4 Natural Sciences    87
5 Social Sciences    116

Two details deserve attention. The order of conditions in case_when() matters: “social sciences” also contains “sciences”, so Social Sciences is checked before Natural Sciences. And across() applies the same change to several columns at once, here turning Yes, yes, and Y into Yes.

After each fix, the categories should be counted again, as above. If a spelling was missed, it shows up immediately.

3.7.5 Missing codes and numbers stored as text

Survey tools often mark missing answers with codes such as 99 or -9. As the opening example showed, these codes must become NA before any calculation, or R will treat them as real values. The function na_if() turns a code into NA, and as.integer() then turns the text into whole numbers:

responses <- responses |>
  mutate(
    financial_worry = financial_worry |> na_if("99") |> na_if("-9") |> as.integer(),
    across(all_of(item_names), ~ as.integer(na_if(.x, "99"))),
    workshop_sessions = as.integer(workshop_sessions),
    age = parse_number(age)
  )

The function parse_number() from readr reads the number out of a text value, ignoring extra text around it. It is useful for the sleep columns, where some students typed "7 hrs". A decimal comma, however, defeats it:

parse_number(c("7 hrs", "6.5", "6,5"))

[1]  7.0  6.5 65.0

The value "6,5" becomes 65, not 6.5, because in English a comma separates thousands. The comma has to be changed to a point first:

parse_number(str_replace(c("7 hrs", "6.5", "6,5"), ",", "."))

[1] 7.0 6.5 6.5

Mistakes like this produce no error message. The only protection is to check the data after each step.

3.7.6 Impossible values

Some values cannot be true, such as an age of 250 or 26 hours of sleep a night. Left in the data, a single one can dominate an average or a regression. They are typing mistakes, and since the true value cannot be known, the honest fix is to set them to NA rather than to guess (an age of 250 might have been 25, or 50):

responses |> filter(age > 100) |> select(student_id, age)

# A tibble: 1 × 2
  student_id   age
  <chr>      <dbl>
1 S0007        250

responses <- responses |> mutate(age = if_else(age > 100, NA, age))

This is where the missing age in Chapter 1 came from. Unusual but possible values, such as a student who sleeps 4 hours a night, are a different matter: they are real data and are kept. Chapter 6 discusses how to handle them.

3.7.7 Scale scores

A questionnaire scale is scored as the average of its items, and the items must all point in the same direction before they are averaged. One stress item, stress_4 (“I feel confident handling problems in my studies”), is worded in the opposite direction to the others: agreeing with it means less stress. Averaged as it stands, it would cancel part of what the other five items measure. It must first be reversed, so that 1 becomes 5, 2 becomes 4, and so on. On a 1-to-5 scale, subtracting from 6 does exactly that:

questionnaire_clean <- responses |>
  select(student_id, all_of(item_names))

scores <- questionnaire_clean |>
  mutate(
    stress_4           = 6 - stress_4,
    stress_score       = rowMeans(pick(stress_1:stress_6), na.rm = TRUE),
    burnout_score      = rowMeans(pick(burnout_1:burnout_6), na.rm = TRUE),
    support_score      = rowMeans(pick(support_1:support_6), na.rm = TRUE),
    satisfaction_score = rowMeans(pick(satisfaction_1:satisfaction_4), na.rm = TRUE)
  ) |>
  select(student_id, ends_with("_score"))
head(scores)

# A tibble: 6 × 5
  student_id stress_score burnout_score support_score satisfaction_score
  <chr>             <dbl>         <dbl>         <dbl>              <dbl>
1 S0191              2.83          1.5           3.83               3.33
2 S0015              4.5           3.83          1.67               3   
3 S0474              3.5           4.5           1.4                2.75
4 S0418              3.17          3             1.67               3.75
5 S0049              3             3.5           3.17               3.33
6 S0537              2.83          2.5           3.2                4.75

The function pick() selects the columns of one scale, and rowMeans() averages each student’s answers across them. With na.rm = TRUE, a student who skipped one item still gets a score from the items they answered. Whether that is acceptable, and how many skipped items is too many, is a decision to report in the thesis. Chapter 9 checks that the items really do measure four separate scales.

The original items should stay in the data, with the scores stored separately as here, because the items are needed again for checking the scales.

3.7.8 From wide to long

The semester measurements sit side by side: GPA_S1, Sleep_S1, …, Wellbeing_S4, 28 columns in all. The data is not tidy, since each row holds four observations, and analyses of change over time need one row per student per semester. This is a pivot_longer() with one extra idea: each column name holds two pieces of information, the variable (GPA) and the semester (1), separated by _S. The special name .value tells R to use the first piece as a column name:

semesters_clean <- responses |>
  select(student_id, matches("_S[1-4]$")) |>
  pivot_longer(
    cols      = -student_id,
    names_to  = c(".value", "semester"),
    names_sep = "_S"
  )
head(semesters_clean)

# A tibble: 6 × 9
  student_id semester GPA   Sleep Study Exercise Caffeine Meetings Wellbeing
  <chr>      <chr>    <chr> <chr> <chr> <chr>    <chr>    <chr>    <chr>    
1 S0191      1        3.07  5.6   31    3        140      2        64       
2 S0191      2        3.29  5,7   30    4        140      7        65       
3 S0191      3        3.35  5.5   32    3        130      5        60       
4 S0191      4        3.58  6.6   29    2        155      7        69       
5 S0015      1        2.8   6.5   15    2        210      2        34       
6 S0015      2        3.2   6.5   18    0        275      5        38       

The selection matches("_S[1-4]$") picks the columns whose names end in _S1 to _S4. The rest is cleaning already familiar from the steps above: tidy names, text to numbers (with the decimal-comma fix for sleep), impossible values to NA, and finally removing the empty rows of students who had left the study:

semesters_clean <- semesters_clean |>
  rename(gpa = GPA, sleep_hours = Sleep, study_hours = Study,
         exercise_days = Exercise, caffeine_mg = Caffeine,
         supervisor_meetings = Meetings, wellbeing = Wellbeing) |>
  mutate(
    semester    = as.integer(semester),
    sleep_hours = parse_number(str_replace(sleep_hours, ",", ".")),
    across(c(gpa, study_hours, exercise_days, caffeine_mg, supervisor_meetings, wellbeing),
           parse_number),
    sleep_hours = if_else(sleep_hours > 24, NA, sleep_hours),
    study_hours = if_else(study_hours > 168, NA, study_hours),
    gpa         = if_else(gpa > 4, NA, gpa)
  ) |>
  filter(!if_all(gpa:wellbeing, is.na))
nrow(semesters_clean)

[1] 2326

The condition if_all(gpa:wellbeing, is.na) is TRUE for rows where every measurement is missing, and ! keeps the others.

3.7.9 Checking the result

The last step checks what the cleaning has produced. The amount of missing data comes first, and across() inside summarise() counts the NAs in every column at once:

semesters_clean |> summarise(across(everything(), ~ sum(is.na(.x))))

# A tibble: 1 × 9
  student_id semester   gpa sleep_hours study_hours exercise_days caffeine_mg
       <int>    <int> <int>       <int>       <int>         <int>       <int>
1          0        0     1          21          23            20          25
# ℹ 2 more variables: supervisor_meetings <int>, wellbeing <int>

Up to 25 values per measurement are still missing, about 1% of the records: answers students really did not give, and the impossible values set to NA. These are genuine gaps, and Chapter 6 discusses how to handle them.

The most satisfying check of all is a comparison with a trusted result. The clean tables in the data2thesis package were produced from this same export, so the cleaned data should match them exactly:

same <- function(mine, theirs) {
  mine   <- as.data.frame(mine)[, names(mine)]
  theirs <- as.data.frame(theirs)[, names(mine)]
  isTRUE(all.equal(mine, theirs, check.attributes = FALSE))
}
same(questionnaire_clean |> arrange(student_id),
     questionnaire |> arrange(student_id))

[1] TRUE

same(semesters_clean |> arrange(student_id, semester),
     semesters |> arrange(student_id, semester))

[1] TRUE

Both are TRUE. The cleaning is complete, and from here on the book uses the clean tables from the package.

write_csv(semesters_clean, here::here("data", "semesters_clean.csv"))

airquality |>
  summarise(
    mean_ozone   = mean(Ozone, na.rm = TRUE),
    missing_days = sum(is.na(Ozone)),
    .by = Month
  )

  Month mean_ozone missing_days
1     5   23.61538            5
2     6   29.44444           21
3     7   59.11538            5
4     8   59.96154            5
5     9   31.44828            1

3.8.1 Summary

• Clean data is valid, accurate, complete as far as possible, consistent, and unique. Cleaning decisions shape the results, so they are made in code from an untouched raw file, recorded, reported, and checked after every step.
• Tidy data has one variable per column, one observation per row, and one value per cell.
• The pipe |> passes a result to the next function: read it as “and then”.
• dplyr verbs: filter() chooses cases, select() chooses variables, arrange() sorts, mutate() creates variables, summarise() with .by summarises by group, and count() counts. case_when() sorts values into categories.
• pivot_longer() and pivot_wider() reshape data between wide and long formats.
• left_join() combines tables by an identifier; anti_join() finds rows without a match.
• Cleaning survey data: remove test and duplicate rows, rename variables, fix categories by pattern, turn missing codes into NA column by column, convert text to numbers (watch decimal commas), set impossible values to NA, reverse reversed items before computing scale scores, and check the result.

3.8.2 Key terms

Data cleaning, tidyverse, tidy data, tibble, pipe, verb, grouped summary, wide format, long format, join, identifier, missing code, reversed item, scale score.

4.5.1 Histograms

The first question about any numeric variable concerns its distribution: where most values lie, how widely they spread, and whether some are unusual. A histogram answers it by dividing the values into ranges, called bins, and showing how many cases fall into each. Here is the distribution of sleep:

R Copy

ggplot(first_sem, aes(x = sleep_hours)) +
  geom_histogram(binwidth = 0.5)

Figure 4.5 shows a roughly symmetrical, bell-shaped distribution, centred a little above 6 hours. Most students sleep less than the recommended 7 hours: 65% of them in the first semester. A few sleep very little, under 4.5 hours; Chapter 6 looks at such unusual values.

The argument binwidth sets the width of each bin, here half an hour. The choice matters: bins that are too wide hide the shape, and bins that are too narrow make it noisy, so it is worth trying a few widths.

4.5.2 Bar charts

A bar chart shows how many cases fall into each category. The geom geom_bar() does the counting:

R Copy

ggplot(students, aes(x = faculty)) +
  geom_bar()

When the heights have already been calculated, for example averages, geom_col() is used instead, with the calculated value mapped to y. Because a bar represents its value by its length, the axis of a bar chart must start at zero; the section on honest graphs returns to this.

4.5.3 Box plots

A box plot summarises a numeric variable for each group, which makes it the standard graph for comparing groups. The box covers the middle half of the values, the line inside is the median, the whiskers reach out to the typical range, and points beyond the whiskers mark unusual values. The comparison below concerns wellbeing in full-time and part-time students:

R Copy

ggplot(first_sem, aes(x = study_mode, y = wellbeing)) +
  geom_boxplot()

Part-time students’ wellbeing is a little lower, but the boxes overlap considerably: many part-time students are doing better than many full-time students. Whether the difference is larger than chance is a question for Chapter 8.

4.5.4 Scatter plots

The relationship between two numeric variables, such as sleep and grades, is shown with a scatter plot, which puts one variable on each axis. With 600 students, points pile on top of each other, a problem called overplotting, so alpha makes them partly transparent, and darker areas show where many students are:

R Copy

ggplot(first_sem, aes(x = sleep_hours, y = gpa)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm")

The layer geom_smooth(method = "lm") adds a straight trend line, with a grey band showing its uncertainty. The line rises: students who sleep more tend to have slightly higher GPAs. The points, however, scatter widely around it. Sleep is related to GPA, but it is far from the whole story, and Chapter 8 measures the relationship.

4.5.5 Line charts

The study followed its students for four semesters, and a line chart shows how the average changed. The averages are calculated first (Chapter 3) and then plotted, with one line per workshop group:

R Copy

wellbeing_by_sem <- semesters |>
  left_join(students, join_by(student_id)) |>
  summarise(mean_wellbeing = mean(wellbeing), .by = c(semester, workshop))

ggplot(wellbeing_by_sem, aes(x = semester, y = mean_wellbeing, colour = workshop)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5)

Figure 4.9 tells the story of the workshop at a glance: the groups start level, the invited students pull ahead in semester 2, and the gap narrows afterwards. Chapters 7 and 10 test whether this pattern is real.

4.7.1 Labels

A reader cannot interpret a graph whose axes say mean_wellbeing and semester. Axis labels should state what was measured, with units, and the title or caption should state the message. The function labs() sets the title, subtitle, axis labels, legend title, and caption:

R Copy

p <- ggplot(wellbeing_by_sem, aes(x = semester, y = mean_wellbeing, colour = workshop)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  labs(
    title    = "Wellbeing over the two years",
    subtitle = "Average score per semester, by workshop group",
    x        = "Semester",
    y        = "Average wellbeing (0 to 100)",
    colour   = "Workshop"
  )
p

Storing the graph in an object, here p, makes it possible to add to it step by step without repeating the code.

4.7.2 Scales and colours

Scales control how data values become positions, colours, and sizes, and each has a scale_ function. Two decisions are needed here. The x axis should show only whole semesters, since there is no semester 2.5. And the colours should come from a palette that people with colour vision deficiency, about one man in twelve, can tell apart. The viridis palettes, built into ggplot2, are designed for exactly that, and they also print well in greyscale:

R Copy

p <- p +
  scale_x_continuous(breaks = 1:4) +
  scale_colour_viridis_d(end = 0.8)
p

Colours can also be chosen by hand, with scale_colour_manual(values = c("Invited" = "#1b9e77", "Not invited" = "#7570b3")). Whatever the choice, colour should never be the only way to tell groups apart in a printed figure: different shapes or line types, or labels placed directly on the lines, keep the graph readable in black and white.

4.7.3 Themes

A theme controls everything that is not data: the background, grid lines, fonts, and the position of the legend. The default grey background of ggplot2 is useful on screen, but most journals prefer a plain look. The complete themes theme_minimal() and theme_bw() provide it, base_size sets the font size, and theme() adjusts individual details:

R Copy

p <- p +
  theme_minimal(base_size = 13) +
  theme(legend.position = "bottom")
p

Compared with Figure 4.9, Figure 4.14 shows the same data, but it is now ready for a thesis.

4.7.4 Annotations

A reference value often helps a reader interpret a graph: a recommended amount, a threshold, a mean. The functions geom_vline() and geom_hline() draw reference lines, and annotate() places text:

R Copy

ggplot(first_sem, aes(x = sleep_hours)) +
  geom_histogram(binwidth = 0.5, fill = "steelblue", colour = "white") +
  geom_vline(xintercept = 7, linetype = "dashed") +
  annotate("text", x = 7.1, y = Inf, label = "Recommended: 7 hours", hjust = 0, vjust = 1.5) +
  labs(x = "Hours of sleep per night", y = "Number of students") +
  theme_minimal(base_size = 13)

In the annotation, y = Inf places the text at the top of the plot, whatever the height of the bars; hjust = 0 makes it start just right of the line, and vjust = 1.5 moves it slightly down from the edge. The line makes the graph’s message, that most students sleep less than recommended, visible without any further explanation.

4.8.1 A graph improved

The principles are clearest when applied to one graph. Suppose the question is whether wellbeing differs between faculties. A first attempt, using defaults and a zoomed axis, might look like Figure 4.16:

R Copy

faculty_means <- first_sem |>
  summarise(wellbeing = mean(wellbeing), .by = faculty)

ggplot(faculty_means, aes(x = faculty, y = wellbeing, fill = faculty)) +
  geom_col() +
  coord_cartesian(ylim = c(58, 62.5))

The graph suggests large differences: the Education bar looks several times taller than the Social Sciences bar. But the axis starts at 58, so the lengths of the bars mean nothing; the averages actually differ by about 3 points on a 100-point scale. The colours add nothing that the axis labels do not already say, the legend repeats them a second time, the faculties appear in alphabetical rather than meaningful order, and the axis titles are variable names. Most seriously, the graph shows only averages, and hides how much students within each faculty differ.

Figure 4.17 shows the same data, redesigned. Every student appears as a faint point, so the spread within each faculty is visible, and the average of each faculty is marked in red. The faculties are ordered by their average, the labels say what was measured, and the legend is gone:

R Copy

ggplot(first_sem, aes(x = wellbeing, y = reorder(faculty, wellbeing))) +
  geom_jitter(height = 0.15, alpha = 0.25, colour = "grey45") +
  stat_summary(fun = mean, geom = "point", size = 3.5, colour = "#b2182b") +
  labs(x = "Wellbeing in semester 1 (0 to 100)", y = NULL) +
  theme_minimal(base_size = 12)

The function reorder() orders the faculties by their average wellbeing, geom_jitter() spreads the points a little vertically so that they do not hide each other, and stat_summary() calculates and draws each faculty’s mean. The redesigned graph tells the truth that the first one hid: the differences between faculties are small compared with the differences between students within each faculty. Chapter 8 tests whether they are larger than chance.

4.9.1 Small multiples

Facets, also called small multiples, split one graph into a panel for each group, all with the same axes, so the groups are easy to compare. The function facet_wrap() takes the variable to split by:

R Copy

ggplot(first_sem, aes(x = sleep_hours, y = gpa)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm") +
  facet_wrap(~ programme) +
  labs(x = "Hours of sleep per night", y = "GPA") +
  theme_minimal(base_size = 12)

The relationship looks similar in both programmes. For two grouping variables, facet_grid(rows ~ columns) makes a grid of panels. Because the panels share their axes, facets respect the rule from the section on reading graphs: comparisons are made by position on a common scale.

4.9.2 Density plots

A density plot is a smoothed histogram. Its advantage is that several distributions can be drawn on top of each other and compared. Here, fill is mapped to the workshop group and alpha lets the two curves show through each other:

R Copy

semesters |>
  filter(semester == 2) |>
  left_join(students, join_by(student_id)) |>
  ggplot(aes(x = wellbeing, fill = workshop)) +
  geom_density(alpha = 0.5) +
  scale_fill_viridis_d(end = 0.8) +
  labs(x = "Wellbeing (0 to 100)", y = "Density", fill = "Workshop") +
  theme_minimal(base_size = 13)

In this code, the data flows straight into ggplot() through the pipe. The two curves have the same shape, but the invited group’s is shifted to the right: the whole distribution moved, not just a few students.

4.9.3 Heat maps

A heat map shows a table of numbers as coloured tiles, and it suits correlation matrices well. Chapter 2 calculated the correlations between four first-semester measurements. To plot them, the matrix is turned into a long table with one row per pair (Chapter 3), and then geom_tile() draws a tile for each pair and geom_text() writes the value on it:

R Copy

library(tidyr)

correlations <- first_sem |>
  select(gpa, sleep_hours, study_hours, wellbeing) |>
  cor(use = "complete.obs")

correlations |>
  as.data.frame() |>
  mutate(var1 = rownames(correlations)) |>
  pivot_longer(-var1, names_to = "var2", values_to = "r") |>
  ggplot(aes(x = var1, y = var2, fill = r)) +
  geom_tile() +
  geom_text(aes(label = round(r, 2))) +
  scale_fill_gradient2(limits = c(-1, 1)) +
  labs(x = NULL, y = NULL, fill = "Correlation") +
  theme_minimal(base_size = 12)

The scale scale_fill_gradient2() is a diverging palette: two colours for negative and positive values, with white at zero, so the direction and strength of each correlation are visible at once. Students who study more hours sleep less (a negative correlation), and students who sleep more report higher wellbeing (positive). The numbers written on the tiles matter, because colour alone is read imprecisely. Chapter 6 explains how to read correlations.

4.12.1 Summary

• Numerical summaries can hide very different patterns, as Anscombe’s datasets show, so data should always be graphed. Exploratory graphs help the researcher understand the data; explanatory graphs show readers a finding.
• People judge positions and lengths most accurately, then angles, areas, and colours. The most important comparison belongs in position or length, within one panel.
• The graph follows from the question and the kinds of variables: histograms for distributions, bar charts for counts, box plots or points for group comparisons, scatter plots for relationships, line charts for change over time.
• Every ggplot2 graph combines data, aesthetic mappings in aes(), and geometric layers, joined with +. Map a colour inside aes() when it should depend on the data; set it outside aes() when it should be fixed.
• labs() adds labels, scale_ functions control axes and colours (viridis palettes are colour-blind-friendly), and themes control the look.
• Honest graphs start bar axes at zero, state when an axis is zoomed, avoid clutter, and use colour appropriately: qualitative palettes for categories, sequential or diverging palettes for quantities.
• Facets split a graph into comparable panels; density plots compare distributions; heat maps show tables of numbers.
• ggsave() exports figures. Set the size in page units, use at least 300 dpi for images, and prefer PDF or SVG when accepted.

4.12.2 Key terms

Exploratory graph, explanatory graph, graphical perception, grammar of graphics, aesthetic mapping, geom, layer, histogram, bin, bar chart, box plot, median, scatter plot, overplotting, trend line, line chart, mapping, setting, scale, theme, annotation, aspect ratio, qualitative palette, sequential palette, diverging palette, facet, density plot, heat map, dpi, vector format.

5.2.1 Three kinds of question

Research questions come in three kinds, and each kind needs a different design and analysis. Descriptive questions ask what is, such as how many hours graduate students sleep; they need a sample that represents the population well. Relational questions ask what goes with what, such as whether students who sleep more have higher GPAs; they need measurements of both variables, and they establish association, not cause. Causal questions ask what leads to what, such as whether the workshop improves wellbeing; they need a design that rules out other explanations, ideally an experiment.

The kind of question decides the kind of claim you can make. Much confusion in published research comes from asking a relational question and answering it in causal words (“sleep improves grades”). The section on study designs below returns to this.

5.3.1 Falsifiability

The philosopher Karl Popper argued that what separates a scientific claim from other claims is not that it can be proved, but that it can be refuted: it forbids certain results (Popper 1959). A claim that is compatible with every possible result tells us nothing. Compare:

A good test of your own hypothesis is to write down, before seeing any data, what result would count against it. If you cannot, the hypothesis is not yet precise enough. If you can, you have also protected yourself against a common temptation: finding, after the fact, a reason why the disappointing result supports your idea after all.

A good hypothesis therefore has four properties. It is about a stated population and stated variables: graduate students at this university, and wellbeing measured with the index. It is testable with the data available, because the variables are measured and there are enough cases. It is falsifiable: it names the result that would contradict it. And it is stated before the data is analysed, since a hypothesis written after seeing the data is a description of the data, not a test of it.

5.3.2 Directional and non-directional hypotheses

A directional hypothesis predicts the direction of an effect (“invited students will have higher wellbeing”). A non-directional hypothesis predicts only that there is a difference (“wellbeing will differ between faculties”). Use a directional hypothesis when theory gives you a clear prediction; use a non-directional one when it does not. Note that the direction of the hypothesis and the choice of a one-tailed or two-tailed test are separate decisions: many researchers state a directional hypothesis but still use a two-tailed test, the cautious choice explained in Chapter 7.

5.3.3 The null hypothesis

Statistical tests do not test the researcher’s hypothesis directly. They test its opposite, the null hypothesis (\(H_0\)): the claim that nothing is going on, no difference and no relationship. The researcher’s hypothesis becomes the alternative hypothesis (\(H_1\)). For the workshop, \(H_1\) states that invited and not-invited students differ in average wellbeing at the end of semester 2, and \(H_0\) that they have the same average wellbeing.

Testing the opposite of what one believes may seem strange, but the null hypothesis has one great advantage: it is precise enough to calculate with. “The workshop helps” does not say how much, but “the workshop makes no difference” makes a sharp prediction: any difference between the groups is due to chance alone, to which students happened to be invited. If we know how large chance differences usually are, we can ask whether the difference in the data is larger than chance would easily produce.

You can see this world of chance by simulation. Suppose the workshop truly does nothing. Then being invited is just a label on students whose wellbeing is unaffected. The code below creates 300 imaginary students with wellbeing scores similar to those in the study (average 60, standard deviation about 11), splits them into two random groups of 150, and records the difference between the group averages. It then repeats this 5,000 times:

set.seed(2026)
chance_differences <- replicate(5000, {
  wellbeing <- rnorm(300, mean = 60, sd = 11)
  group     <- sample(rep(c("Invited", "Not invited"), 150))
  mean(wellbeing[group == "Invited"]) - mean(wellbeing[group == "Not invited"])
})
summary(chance_differences)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-5.188460 -0.871667  0.001230 -0.008296  0.838404  4.688888 

The function replicate() runs the code in curly brackets many times and collects the results, rnorm() draws random numbers from a normal distribution, and sample() shuffles the labels. Figure 5.2 shows the 5,000 differences:

R Copy

library(ggplot2)
ggplot(data.frame(difference = chance_differences), aes(x = difference)) +
  geom_histogram(binwidth = 0.5, fill = "grey70", colour = "white") +
  geom_vline(xintercept = quantile(chance_differences, c(0.025, 0.975)),
             linetype = "dashed") +
  labs(x = "Difference in average wellbeing (invited minus not invited)",
       y = "Number of simulated studies")

Even when the workshop does nothing, the two groups almost never have exactly the same average: chance alone produces differences, usually small ones. In 95% of the simulated studies, the difference lies between -2.5 and 2.6 points (the dashed lines). A real study that found a difference of 5 points would therefore be hard to explain by chance, and would count as evidence against \(H_0\). A difference of 1 point would not: it is exactly what chance produces. Chapter 7 turns this idea into the p-value, using the real data.

Two things follow from this picture. First, rejecting \(H_0\) is not proving \(H_1\): it says only that “nothing is going on” is a poor explanation of the data. Second, not rejecting \(H_0\) is not proving it: a small study can fail to detect a real effect, just as a blurred photograph can fail to show a real face.

5.3.4 Hypotheses for the whole study

Table 5.4 writes out the study’s research questions in this form. Some questions are confirmatory: they test a hypothesis stated in advance. Others are exploratory: they look for patterns without a prior prediction, and their results suggest hypotheses for future studies rather than confirming them. Both are legitimate, as long as they are labelled honestly.

The later chapters return to this table: each restates its hypothesis before the analysis, and ends by saying whether the data contradicts the null hypothesis.

5.4.1 Cases and the unit of analysis

As Chapter 2 showed, data is a table of cases (rows) and variables (columns). The first question about any dataset is what one row represents: the unit of analysis. The wellbeing study has three, in different tables:

nrow(students)      # one row per student

[1] 600

nrow(semesters)     # one row per student per semester

[1] 2326

nrow(supervisors)   # one row per supervisor

[1] 120

The unit of analysis decides what a question means. “Do students who sleep more have higher GPAs?” is a question about students, so each student should count once, for example with their first-semester values or their average over the semesters. Treating the 2326 semester rows as 2326 independent cases would count each student up to four times, and make the evidence look stronger than it is. Chapter 10 shows how to use all the rows correctly.

5.4.2 From constructs to variables

Many things researchers care about cannot be observed directly: stress, wellbeing, motivation, intelligence, job satisfaction, quality of life. Such ideas are called constructs. To study a construct, you must decide how to measure it, a step called operationalisation. Every operationalisation is a choice, and a different choice could give a different answer.

Table 5.5 shows how the wellbeing study operationalises some of its constructs.

Stress is a good example. No single question captures it, so the questionnaire asks six, each about one aspect, and averages the answers. Look at the items, and at how the answers of the first few students differ from item to item:

library(dplyr)
questionnaire |>
  select(student_id, stress_1:stress_6) |>
  head(5)

  student_id stress_1 stress_2 stress_3 stress_4 stress_5 stress_6
1      S0001        3        2        4        4        3        5
2      S0002        4        3       NA        2        5        3
3      S0003        4        4        2        1        4        3
4      S0004        4        5        5        1        4        5
5      S0005        4        4        4        3        4        2

One item, stress_4 (“I feel confident handling problems in my studies”), is worded in the opposite direction: agreeing with it means less stress. Before averaging, it must be reversed, so that 5 becomes 1 and 1 becomes 5. You did this in Chapter 3:

stress <- questionnaire |>
  mutate(stress_4 = 6 - stress_4,
         stress   = rowMeans(pick(stress_1:stress_6), na.rm = TRUE)) |>
  select(student_id, stress)
summary(stress$stress)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.333   2.667   3.167   3.216   3.667   5.000 

The resulting score is the study’s operational definition of stress. Anyone who reads the thesis should be able to see exactly how it was built, which is why the methods chapter of a thesis describes every measure, and why the questionnaire items go in an appendix.

5.4.3 Levels of measurement

Not all numbers are numbers in the same sense. The psychologist S. S. Stevens proposed four levels of measurement, each allowing more than the one before (Stevens 1946):

The level decides which summaries and tests make sense. The average faculty is meaningless. “Twice as much caffeine” is meaningful (caffeine has a true zero); “twice as much wellbeing” is not (a score of 0 does not mean no wellbeing at all).

R does not know the level of measurement; you have to tell it. Look at how R stores employment:

table(students$employment)

Full-time job          None Part-time job 
           97           307           196 

R lists the categories in alphabetical order, which puts “Full-time job” before “None”. For a nominal variable, that would not matter. But employment is ordinal: none, part-time, full-time is a real order, and tables and graphs should show it. Tell R with an ordered factor:

students <- students |>
  mutate(employment = factor(employment,
                             levels = c("None", "Part-time job", "Full-time job"),
                             ordered = TRUE))
table(students$employment)

         None Part-time job Full-time job 
          307           196            97 

Now tables and graphs follow the natural order, and comparisons such as employment > "None" work.

The opposite mistake is more common: an ordinal variable stored as numbers, such as financial_worry (1 = not at all worried, 5 = extremely worried). R will happily calculate its mean, but the distance from “not at all” to “a little” need not equal the distance from “very” to “extremely”. For a single ordinal item, report the distribution or the median:

table(students$financial_worry, useNA = "ifany")

   1    2    3    4    5 <NA> 
  84  139  172  111   56   38 

median(students$financial_worry, na.rm = TRUE)

[1] 3

Scale scores, the average of several items, are usually treated as interval data, and the book does so. Averaging several items smooths out the unequal steps of each one. This is a convention with good support, not a law; if an examiner asks, say that you followed it, and why.

5.4.4 The roles variables play

In an analysis, variables also have roles. The outcome, or dependent variable, is what the researcher wants to explain or predict, such as wellbeing, GPA, or considering dropout. A predictor, also called an independent or explanatory variable, is what is used to explain or predict it, such as the workshop, sleep, or stress. Three further roles concern a third variable that affects the relationship between the two. A confounder is related to both the predictor and the outcome, and can create a misleading association between them. A moderator changes the strength or direction of a relationship, so that the effect of the predictor depends on it. A mediator lies on the path between predictor and outcome: it is how the predictor has its effect.

The roles are not properties of the variables; they come from the research question. Stress is an outcome in one question (whether the workshop reduces stress) and a predictor in another (whether stress predicts dropout).

Figure 5.3 draws the three less familiar roles as diagrams, using examples from the study: arrows show which variable is thought to influence which.

flowchart LR
  subgraph Confounder
    S1[Sleep] --> C1[Caffeine]
    S1 --> G1[GPA]
    C1 -. "apparent link" .- G1
  end
  subgraph Moderator
    SU[Support] --> G2[GPA]
    P[Programme:<br/>Master's or PhD] --> M((" "))
    M -.-> G2
  end
  subgraph Mediator
    W[Workshop] --> ST[Lower stress] --> WB[Wellbeing]
  end

The confounder is the most important of the three, because it can produce an association that is not a cause. In the wellbeing data, students who take more caffeine have lower GPAs. But students who take more caffeine also sleep less: