Welcome

This notebook is the hands-on companion to the 5-day R & RStudio for Health Sciences course. It is self-contained: the setup chunk below generates a small synthetic clinical-trial dataset so every chunk in this notebook runs top to bottom with no external files required. In a real course, data-raw/ would instead hold your own de-identified export (REDCap, Excel, SPSS, Stata, …) — everything you learn here about importing and cleaning transfers directly.

How to use this notebook in RStudio: open the .Rmd file, then run chunks one at a time with the green ▶ button (or Ctrl+Shift+Enter), or click Preview to render the whole thing to .nb.html. Text like this paragraph is Markdown; grey boxes are R code chunks; output appears directly underneath each chunk.

Project setup

# Packages used across the week. In a normal environment you would run
# install.packages() once for anything missing, then library() every
# session. This chunk loads only the packages that ship with a standard
# tidyverse + health-stats setup.
library(tidyverse)   # dplyr, tidyr, ggplot2, readr, stringr, purrr, forcats
library(janitor)      # clean_names() and friends
library(lubridate)     # dates
library(survival)       # Surv(), survfit(), coxph() -- ships with base R
library(broom)            # tidy() -- turns model objects into data frames
library(haven)              # SPSS / Stata import
library(readxl)              # Excel import
library(flextable)            # Word-ready tables

knitr::opts_chunk$set(fig.width = 7, fig.height = 4.5, dpi = 150)

dir.create("data-raw",   showWarnings = FALSE)
dir.create("data-clean", showWarnings = FALSE)
dir.create("outputs",    showWarnings = FALSE)

Generating the practice dataset (run once)

This chunk plays the role your own 00_generate_sample_data.R would play on a real project: it creates three CSV files inside data-raw/, exactly as if they had been exported from REDCap / Excel. Everything from Day 1 onward reads these files back in, the same way you would read your own data.

set.seed(2026)
n <- 220

sites <- c("Site A", "Site B", "Site C")
arms  <- c("Control", "Intervention")

demographics <- data.frame(
  patient_id   = sprintf("P%03d", 1:n),
  age          = pmin(pmax(round(rnorm(n, 52, 14)), 18), 90),
  sex          = sample(c("F", "M"), n, replace = TRUE, prob = c(0.55, 0.45)),
  site         = sample(sites, n, replace = TRUE),
  treatment_arm = sample(arms, n, replace = TRUE),
  weight_kg    = round(rnorm(n, 75, 14), 1),
  height_m     = round(rnorm(n, 1.68, 0.09), 2),
  icd10        = sample(c("I10","I25","E11","E11.9","J45","J44","M54","N39"),
                         n, replace = TRUE),
  satisfaction_raw = sample(
    c("Very dissatisfied","Dissatisfied","Neutral","Satisfied","Very satisfied"),
    n, replace = TRUE, prob = c(0.05,0.10,0.20,0.40,0.25)),
  time_months  = pmin(round(rexp(n, rate = 1/18), 1), 36),
  stringsAsFactors = FALSE
)

demographics$event <- rbinom(n, 1, ifelse(demographics$treatment_arm == "Intervention", 0.28, 0.42))

bmi_tmp <- demographics$weight_kg / (demographics$height_m^2)
p_diab  <- plogis(-6 + 0.04 * demographics$age + 0.09 * bmi_tmp)
demographics$diabetic <- rbinom(n, 1, p_diab) == 1

demographics$outcome <- round(rnorm(n, 5, 2), 1)
demographics$outcome[sample(seq_len(n), round(0.06 * n))] <- NA   # a little missingness, on purpose

messy_idx <- sample(seq_len(n), round(0.08 * n))                    # a little messiness, on purpose
demographics$site[messy_idx] <- paste0("  ", tolower(demographics$site[messy_idx]), " ")

write.csv(demographics, "data-raw/demographics.csv", row.names = FALSE)

labs1 <- data.frame(
  patient_id = demographics$patient_id, visit = 1,
  sbp = round(rnorm(n, 128 + 8 * demographics$diabetic, 12)),
  dbp = round(rnorm(n, 80, 8)))
dup_ids <- sample(demographics$patient_id, round(0.15 * n))          # a few repeat visits, on purpose
labs2 <- data.frame(
  patient_id = dup_ids, visit = 2,
  sbp = round(rnorm(length(dup_ids), 124, 12)),
  dbp = round(rnorm(length(dup_ids), 78, 8)))
write.csv(rbind(labs1, labs2), "data-raw/lab_results.csv", row.names = FALSE)

trial_raw <- data.frame(
  `Patient ID` = demographics$patient_id, `Age` = demographics$age,
  `Sex` = ifelse(demographics$sex == "F", "Female", "Male"),
  `SBP (mmHg)` = labs1$sbp, `Diabetic` = demographics$diabetic,
  check.names = FALSE)
trial_raw[3, "SBP (mmHg)"] <- 0     # an implausible value, on purpose
trial_raw[7, "Age"] <- 199            # an implausible value, on purpose
write.csv(trial_raw, "data-raw/clinical_trial_raw.csv", row.names = FALSE)

cat("data-raw/ ready:", nrow(demographics), "patients,",
    nrow(labs1) + nrow(labs2), "lab visits.\n")
## data-raw/ ready: 220 patients, 253 lab visits.

Day 1 — Foundations (Mon, Sep 21)

Day 1 builds every core R skill from zero: the calculator, the four data structures you’ll use constantly (vectors, lists, matrices, data frames), factors, and the programming basics (if/else, loops, functions) that make everything later possible.

R as a calculator

2 + 2          # addition
## [1] 4
10 - 3         # subtraction
## [1] 7
6 * 7          # multiplication
## [1] 42
100 / 8        # division
## [1] 12.5
2 ^ 5          # exponent -> 32
## [1] 32
17 %% 5        # modulo (remainder) -> 2
## [1] 2
17 %/% 5       # integer division -> 3
## [1] 3
round(3.14159, 2)
## [1] 3.14

Hello, World!

print(), cat() and message() all display text, but they are not interchangeable: print() shows an R object (with quotes), cat() prints raw formatted text, and message() writes to a separate status stream.

print("Hello, World!")
## [1] "Hello, World!"
cat("Hello, World!\n")
## Hello, World!
name <- "Trainee"
paste("Hello,", name)          # paste() glues pieces together, space-separated
## [1] "Hello, Trainee"
paste0("Hello, ", name, "!")    # paste0() = no automatic space
## [1] "Hello, Trainee!"

Order of operations & comparisons

2 + 3 * 4        # 14 -- multiplication before addition
## [1] 14
(2 + 3) * 4      # 20 -- parentheses override
## [1] 20
-2^2             # -4 -- unary minus applies AFTER ^ (a classic gotcha!)
## [1] -4
age <- 45
age >= 65        # is this patient a senior? -> FALSE
## [1] FALSE
5 == 5           # equal to (do NOT confuse with the single = of assignment)
## [1] TRUE

Logical operators

&/| are vectorized (check every element — use inside filter()); &&/|| check a single value and short-circuit (use inside if()).

TRUE & FALSE         # AND (vectorized) -> FALSE
## [1] FALSE
TRUE | FALSE          # OR  (vectorized) -> TRUE
## [1] TRUE
!TRUE                  # NOT              -> FALSE
## [1] FALSE
xor(TRUE, FALSE)         # exclusive OR    -> TRUE
## [1] TRUE
sbp_val <- 145
sbp_val > 140 && sbp_val < 180   # scalar AND, for if()
## [1] TRUE

Data types: typeof() vs class()

typeof(2.5)             # "double"
## [1] "double"
typeof(2L)                # "integer"
## [1] "integer"
typeof("hello")             # "character"
## [1] "character"
typeof(TRUE)                  # "logical"
## [1] "logical"
class(as.Date("2024-01-01"))    # "Date" -- class() is usually what you want
## [1] "Date"

Special values & type coercion

NA = missing, NULL = nothing at all, NaN = undefined math, Inf = infinity. R will silently convert types to the most flexible one in a mix — always confirm with class() after importing real data.

c(1, "2", 3)            # all become character
## [1] "1" "2" "3"
c(TRUE, 1, 2)             # all become numeric (TRUE -> 1)
## [1] 1 1 2
raw_vals <- c("120", "135", "missing", "142")
as.numeric(raw_vals)       # 120 135 NA 142 -- with a warning! Always check.
## Warning: NAs introduced by coercion
## [1] 120 135  NA 142

Vectors: R’s basic building block

sbp <- c(128, 135, 142, 119, 151, 133)
sbp[3]              # 142 -- third patient (R indexing starts at 1!)
## [1] 142
sbp[sbp > 130]      # logical indexing -- only readings above 130
## [1] 135 142 151 133
mean(sbp); sd(sbp); length(sbp)
## [1] 134.6667
## [1] 11.0755
## [1] 6
sbp + 10             # vectorized arithmetic -- no loop needed
## [1] 138 145 152 129 161 143

Sequences, lists & the [ ] vs [[ ]] rule

1:10                          # integers 1 through 10
##  [1]  1  2  3  4  5  6  7  8  9 10
seq(0, 100, by = 25)           # 0 25 50 75 100
## [1]   0  25  50  75 100
patient1 <- list(id = "P001", age = 54, diagnoses = c("HTN", "T2DM"))
patient1$age            # -> 54          (the value itself)
## [1] 54
patient1["age"]           # -> a LIST containing age  (still wrapped)
## $age
## [1] 54
patient1[["age"]]           # -> 54          (same as $, unwrapped)
## [1] 54

Matrices & data frames

lab_matrix <- matrix(c(120, 80, 135, 85, 142, 90), nrow = 3,
  dimnames = list(c("P1","P2","P3"), c("SBP","DBP")))
lab_matrix["P2", "SBP"]     # index like [row, col] -> 135
## [1] 80
patients <- data.frame(id = c("P01","P02","P03"), age = c(54,67,45),
                        sex = c("F","M","F"), sbp = c(128,142,119))
str(patients)
## 'data.frame':    3 obs. of  4 variables:
##  $ id : chr  "P01" "P02" "P03"
##  $ age: num  54 67 45
##  $ sex: chr  "F" "M" "F"
##  $ sbp: num  128 142 119
patients[patients$age > 50, ]     # every patient over 50

Factors: categorical data done properly

severity <- factor(c("mild","severe","moderate","mild"),
  levels = c("mild","moderate","severe"), ordered = TRUE)
levels(severity)
## [1] "mild"     "moderate" "severe"
severity[1] < severity[2]         # TRUE -- ordered comparison!
## [1] TRUE
table(severity)
## severity
##     mild moderate   severe 
##        2        1        1

Control flow, loops & functions

bmi_val <- 31
if (bmi_val < 18.5) {
  "Underweight"
} else if (bmi_val < 25) {
  "Normal weight"
} else if (bmi_val < 30) {
  "Overweight"
} else {
  "Obese"
}
## [1] "Obese"
for (i in 1:5) cat("Square of", i, "is", i^2, "\n")
## Square of 1 is 1 
## Square of 2 is 4 
## Square of 3 is 9 
## Square of 4 is 16 
## Square of 5 is 25
bmi_calc <- function(weight_kg, height_m) weight_kg / (height_m^2)
bmi_calc(70, 1.75)               # 22.86
## [1] 22.85714

Missing data

mean(c(120, 135, NA))                 # NA -- contagious!
## [1] NA
mean(c(120, 135, NA), na.rm = TRUE)   # 127.5
## [1] 127.5

LAB 1: Import, inspect & clean

trial_raw <- read_csv("data-raw/clinical_trial_raw.csv", show_col_types = FALSE)

trial_day1 <- trial_raw |> clean_names()   # "Patient ID" -> patient_id, "SBP (mmHg)" -> sbp_mm_hg
names(trial_day1)
## [1] "patient_id" "age"        "sex"        "sbp_mm_hg"  "diabetic"
dim(trial_day1)
## [1] 220   5
summary(trial_day1)
##   patient_id             age             sex              sbp_mm_hg    
##  Length:220         Min.   : 18.00   Length:220         Min.   :  0.0  
##  Class :character   1st Qu.: 42.75   Class :character   1st Qu.:121.0  
##  Mode  :character   Median : 52.00   Mode  :character   Median :128.0  
##                     Mean   : 53.15                      Mean   :128.0  
##                     3rd Qu.: 61.25                      3rd Qu.:135.2  
##                     Max.   :199.00                      Max.   :158.0  
##   diabetic      
##  Mode :logical  
##  FALSE:174      
##  TRUE :46       
##                 
##                 
## 
saveRDS(trial_day1, "data-clean/trial_day1.rds")

Checkpoint: How many patients are in the dataset? Which variable has implausible values (look closely at age and sbp_mm_hg above — one row of each was deliberately corrupted for this exercise)?

Day 2 — Data Wrangling for Research (Tue, Sep 22)

The pipe

trial_day1 <- readRDS("data-clean/trial_day1.rds")

trial_day1 |>
  filter(age >= 18) |>
  select(age, sbp_mm_hg) |>
  summary()
##       age           sbp_mm_hg    
##  Min.   : 18.00   Min.   :  0.0  
##  1st Qu.: 42.75   1st Qu.:121.0  
##  Median : 52.00   Median :128.0  
##  Mean   : 53.15   Mean   :128.0  
##  3rd Qu.: 61.25   3rd Qu.:135.2  
##  Max.   :199.00   Max.   :158.0

filter(), select(), arrange()

trial_day1 |> filter(age >= 18) |> nrow()                 # adults only
## [1] 220
trial_day1 |> select(patient_id, age, sex) |> head(3)      # keep some columns
trial_day1 |> arrange(desc(age)) |> head(3)                  # oldest first

Preparing for a real join

demographics <- read_csv("data-raw/demographics.csv", show_col_types = FALSE) |> clean_names()
labs         <- read_csv("data-raw/lab_results.csv",  show_col_types = FALSE) |> clean_names()

class(demographics$patient_id); class(labs$patient_id)     # 1. types match?
## [1] "character"
## [1] "character"
labs |> count(patient_id) |> filter(n > 1) |> nrow()          # 2. duplicate IDs?
## [1] 33

Resolving duplicate IDs before joining

labs_latest <- labs |>
  group_by(patient_id) |>
  slice_max(visit, n = 1, with_ties = FALSE) |>
  ungroup()

nrow(labs); nrow(labs_latest)   # confirm duplicates are gone
## [1] 253
## [1] 220

left_join() and mutate(): deriving BMI

combined <- demographics |> left_join(labs_latest, by = "patient_id")
nrow(demographics); nrow(combined)   # 3. did the row count change unexpectedly?
## [1] 220
## [1] 220
combined <- combined |> mutate(bmi = weight_kg / (height_m^2))
combined |> select(patient_id, weight_kg, height_m, bmi) |> head(3)

Age groups with cut()

combined <- combined |>
  mutate(age_group = cut(age, breaks = c(0, 18, 40, 65, Inf),
           labels = c("<18","18-39","40-64","65+"), right = FALSE))
table(combined$age_group)
## 
##   <18 18-39 40-64   65+ 
##     0    39   140    41

group_by() + summarise()

Notice the messy site values (" site a ", "site b", …) are still splitting patients into extra, spurious groups here — this is exactly the problem the string-cleaning step fixes further down.

combined |>
  group_by(site, sex) |>
  summarise(n = n(), mean_age = mean(age, na.rm = TRUE),
            mean_sbp = mean(sbp, na.rm = TRUE), .groups = "drop")

Detecting & handling missing data

colSums(is.na(combined))
##       patient_id              age              sex             site 
##                0                0                0                0 
##    treatment_arm        weight_kg         height_m            icd10 
##                0                0                0                0 
## satisfaction_raw      time_months            event         diabetic 
##                0                0                0                0 
##          outcome            visit              sbp              dbp 
##               13                0                0                0 
##              bmi        age_group 
##                0                0
n_before <- nrow(combined)
combined <- combined |> drop_na(outcome)
n_before - nrow(combined)   # patients excluded -- report this number in your methods
## [1] 13

Recoding with case_when()

combined <- combined |>
  mutate(bp_category = case_when(
    sbp < 120              ~ "Normal",
    sbp >= 120 & sbp < 130 ~ "Elevated",
    sbp >= 130 & sbp < 140 ~ "Stage 1",
    sbp >= 140             ~ "Stage 2",
    TRUE                   ~ NA_character_))
table(combined$bp_category, useNA = "ifany")
## 
## Elevated   Normal  Stage 1  Stage 2 
##       69       44       60       34
combined <- combined |>
  mutate(diagnosis_group = case_when(
    str_starts(icd10, "I")  ~ "Cardiovascular",
    str_starts(icd10, "E1") ~ "Diabetes",
    str_starts(icd10, "J")  ~ "Respiratory",
    TRUE                    ~ "Other"))
table(combined$diagnosis_group)
## 
## Cardiovascular       Diabetes          Other    Respiratory 
##             49             58             50             50

Ordered factors: Likert scales

combined <- combined |>
  mutate(satisfaction = factor(satisfaction_raw,
    levels = c("Very dissatisfied","Dissatisfied","Neutral","Satisfied","Very satisfied"),
    ordered = TRUE))
summary(combined$satisfaction)
## Very dissatisfied      Dissatisfied           Neutral         Satisfied 
##                 8                25                46                69 
##    Very satisfied 
##                59

Bonus: string cleaning with stringr

messy_site <- c("  Site A ", "site b", "SITE C", "Site_A")
str_to_title(str_trim(messy_site))     # consistent Title Case
## [1] "Site A" "Site B" "Site C" "Site_a"
combined <- combined |> mutate(site = str_to_upper(str_trim(site)))
table(combined$site)     # now exactly 3 clean groups
## 
## SITE A SITE B SITE C 
##     66     75     66

Bonus: reshaping with pivot_longer()

wide_example <- data.frame(patient_id = c("P01","P02"),
  sbp_visit1 = c(128, 142), sbp_visit2 = c(122, 138))

wide_example |>
  pivot_longer(cols = starts_with("sbp_visit"), names_to = "visit", values_to = "sbp")

LAB 2 checkpoint & save

trial <- combined
saveRDS(trial, "data-clean/trial_day2.rds")
dim(trial)
## [1] 207  21

Checkpoint: How many patients were excluded for missing outcome? Which diagnosis_group is most common? Compare the group_by(site) output before and after the string-cleaning step — how many groups collapsed into one?

Day 3 — Visualisation & Descriptive Statistics (Wed, Sep 23)

Anatomy of a ggplot() call

Every ggplot needs three ingredients: data, an aes() mapping (which columns go to x/y/colour), and a geom_*() layer for how to draw it.

trial <- readRDS("data-clean/trial_day2.rds")

ggplot(trial, aes(x = age, y = sbp)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(title = "Systolic BP by Age", x = "Age (years)", y = "SBP (mmHg)") +
  theme_minimal()

A consistent, publication-ready theme

my_theme <- theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank(),
        legend.position = "bottom")

Histogram & boxplot

fig1 <- ggplot(trial, aes(x = sbp)) +
  geom_histogram(binwidth = 5, fill = "#16A3A0", colour = "white") +
  labs(title = "Distribution of Systolic Blood Pressure") + my_theme
fig1

fig2 <- ggplot(trial, aes(treatment_arm, sbp, fill = treatment_arm)) +
  geom_boxplot() +
  labs(title = "SBP by Treatment Arm", x = NULL) + my_theme +
  theme(legend.position = "none")
fig2

Bar chart & faceted scatterplot

ggplot(trial, aes(diagnosis_group)) +
  geom_bar(fill = "#FF6B4A") + coord_flip() +
  labs(title = "Patients by Diagnosis Group", x = NULL, y = "Count") + my_theme

ggplot(trial, aes(age, sbp, colour = sex)) +
  geom_point(alpha = 0.6) +
  facet_wrap(~ site) +
  labs(title = "Age vs. SBP by Study Site") + my_theme

Kaplan-Meier survival curves

We build the curve with base survival, then hand it to broom::tidy() so it plots with the same ggplot() grammar as every other figure today — no extra plotting package required.

fit <- survfit(Surv(time_months, event) ~ treatment_arm, data = trial)

km_df <- broom::tidy(fit) |>
  mutate(strata = str_remove(strata, "treatment_arm="))

ggplot(km_df, aes(x = time, y = estimate, colour = strata, fill = strata)) +
  geom_step(linewidth = 1) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.15, colour = NA) +
  labs(x = "Months", y = "Survival probability", colour = "Arm", fill = "Arm",
       title = "Kaplan-Meier Survival by Treatment Arm") +
  my_theme

Table 1: baseline characteristics

gtsummary::tbl_summary() is the standard tool for this if it’s installed (install.packages("gtsummary")); this chunk uses it when available and otherwise falls back to a plain dplyr summary so the notebook always runs.

if (requireNamespace("gtsummary", quietly = TRUE)) {
  library(gtsummary)
  trial |>
    select(age, sex, bmi, sbp, diagnosis_group, treatment_arm) |>
    tbl_summary(by = treatment_arm) |>
    add_p() |> add_overall() |> bold_labels()
} else {
  trial |>
    group_by(treatment_arm) |>
    summarise(n = n(),
              mean_age = round(mean(age, na.rm = TRUE), 1),
              mean_bmi = round(mean(bmi, na.rm = TRUE), 1),
              mean_sbp = round(mean(sbp, na.rm = TRUE), 1),
              .groups = "drop")
}
## 
## Attaching package: 'gtsummary'
## The following object is masked _by_ '.GlobalEnv':
## 
##     trial
## The following object is masked from 'package:flextable':
## 
##     continuous_summary
Characteristic Overall
N = 2071		 Control
N = 1081		 Intervention
N = 991		 p-value2
age 51 (42, 61) 51 (42, 62) 52 (42, 60) 0.9
sex


0.8
    F 113 (55%) 58 (54%) 55 (56%)
    M 94 (45%) 50 (46%) 44 (44%)
bmi 25.9 (22.9, 30.2) 26.6 (22.7, 30.9) 25.8 (23.3, 29.8) 0.5
sbp 128 (121, 137) 128 (121, 135) 128 (121, 138) 0.6
diagnosis_group


0.3
    Cardiovascular 49 (24%) 31 (29%) 18 (18%)
    Diabetes 58 (28%) 26 (24%) 32 (32%)
    Other 50 (24%) 24 (22%) 26 (26%)
    Respiratory 50 (24%) 27 (25%) 23 (23%)
1 Median (Q1, Q3); n (%)
2 Wilcoxon rank sum test; Pearson’s Chi-squared test

Exporting figures & tables

ggsave("outputs/fig1_sbp_distribution.png", fig1, width = 7, height = 5, dpi = 300)
ggsave("outputs/fig2_sbp_by_arm.png", fig2, width = 7, height = 5, dpi = 300)

Checkpoint: Is SBP roughly normal or skewed (fig1)? Which arm has the higher median SBP (fig2)? Which baseline variable, if any, differs between arms in Table 1?

Day 4 — Statistical Analysis for Health Research (Thu, Sep 24)

t-test, chi-square & ANOVA

t.test(sbp ~ treatment_arm, data = trial)                        # 2 groups, continuous
## 
##  Welch Two Sample t-test
## 
## data:  sbp by treatment_arm
## t = -0.89284, df = 204.47, p-value = 0.373
## alternative hypothesis: true difference in means between group Control and group Intervention is not equal to 0
## 95 percent confidence interval:
##  -4.563969  1.718851
## sample estimates:
##      mean in group Control mean in group Intervention 
##                   127.5370                   128.9596
chisq.test(table(trial$diagnosis_group, trial$sex))               # 2 categorical variables
## 
##  Pearson's Chi-squared test
## 
## data:  table(trial$diagnosis_group, trial$sex)
## X-squared = 2.1081, df = 3, p-value = 0.5503
fit_aov <- aov(sbp ~ site, data = trial)                             # 3+ group means
summary(fit_aov)
##              Df Sum Sq Mean Sq F value Pr(>F)
## site          2    169   84.43    0.64  0.528
## Residuals   204  26902  131.87

Non-parametric alternatives

wilcox.test(sbp ~ treatment_arm, data = trial)   # instead of t.test, if skewed
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  sbp by treatment_arm
## W = 5099, p-value = 0.5667
## alternative hypothesis: true location shift is not equal to 0
kruskal.test(sbp ~ site, data = trial)             # instead of ANOVA, if skewed
## 
##  Kruskal-Wallis rank sum test
## 
## data:  sbp by site
## Kruskal-Wallis chi-squared = 2.1094, df = 2, p-value = 0.3483

Linear regression

fit_lm <- lm(sbp ~ age + bmi, data = trial)
summary(fit_lm)
## 
## Call:
## lm(formula = sbp ~ age + bmi, data = trial)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.442  -7.033  -0.009   8.361  31.786 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 124.95360    4.75851  26.259   <2e-16 ***
## age           0.09116    0.05812   1.569    0.118    
## bmi          -0.05553    0.13120  -0.423    0.673    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.45 on 204 degrees of freedom
## Multiple R-squared:  0.01291,    Adjusted R-squared:  0.00323 
## F-statistic: 1.334 on 2 and 204 DF,  p-value: 0.2658
confint(fit_lm)
##                    2.5 %      97.5 %
## (Intercept) 115.57142068 134.3357749
## age          -0.02342837   0.2057386
## bmi          -0.31422049   0.2031561

Logistic regression: odds ratios

fit_glm <- glm(diabetic ~ age + bmi + sex, data = trial, family = binomial)
summary(fit_glm)
## 
## Call:
## glm(formula = diabetic ~ age + bmi + sex, family = binomial, 
##     data = trial)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -6.76685    1.25402  -5.396 6.81e-08 ***
## age          0.05980    0.01445   4.138 3.50e-05 ***
## bmi          0.07535    0.02889   2.608  0.00911 ** 
## sexM         0.05833    0.37317   0.156  0.87579    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 208.82  on 206  degrees of freedom
## Residual deviance: 182.99  on 203  degrees of freedom
## AIC: 190.99
## 
## Number of Fisher Scoring iterations: 5
exp(cbind(OR = coef(fit_glm), confint(fit_glm)))   # exponentiate: log-odds -> odds ratio
## Waiting for profiling to be done...
##                      OR        2.5 %     97.5 %
## (Intercept) 0.001151316 8.571301e-05 0.01201421
## age         1.061623305 1.032990e+00 1.09355291
## bmi         1.078266616 1.019795e+00 1.14299198
## sexM        1.060062446 5.070877e-01 2.20719644

Log-rank test & Cox regression

survdiff(Surv(time_months, event) ~ treatment_arm, data = trial)
## Call:
## survdiff(formula = Surv(time_months, event) ~ treatment_arm, 
##     data = trial)
## 
##                              N Observed Expected (O-E)^2/E (O-E)^2/V
## treatment_arm=Control      108       47     44.2     0.182     0.501
## treatment_arm=Intervention  99       25     27.8     0.289     0.501
## 
##  Chisq= 0.5  on 1 degrees of freedom, p= 0.5
cox_fit <- coxph(Surv(time_months, event) ~ treatment_arm + age + sex, data = trial)
summary(cox_fit)
## Call:
## coxph(formula = Surv(time_months, event) ~ treatment_arm + age + 
##     sex, data = trial)
## 
##   n= 207, number of events= 72 
## 
##                                coef exp(coef)  se(coef)      z Pr(>|z|)  
## treatment_armIntervention -0.158934  0.853053  0.250891 -0.633   0.5264  
## age                       -0.022691  0.977565  0.008816 -2.574   0.0101 *
## sexM                       0.153167  1.165519  0.237340  0.645   0.5187  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                           exp(coef) exp(-coef) lower .95 upper .95
## treatment_armIntervention    0.8531      1.172    0.5217    1.3949
## age                          0.9776      1.023    0.9608    0.9946
## sexM                         1.1655      0.858    0.7320    1.8559
## 
## Concordance= 0.607  (se = 0.042 )
## Likelihood ratio test= 7.79  on 3 df,   p=0.05
## Wald test            = 7.76  on 3 df,   p=0.05
## Score (logrank) test = 7.79  on 3 df,   p=0.05
cox.zph(cox_fit)   # check the proportional-hazards assumption
##                chisq df    p
## treatment_arm 0.0152  1 0.90
## age           1.3124  1 0.25
## sex           0.2358  1 0.63
## GLOBAL        1.4144  3 0.70

Checkpoint: Was the t-test significant? Is BMI a significant predictor of SBP in fit_lm? Which variable has the largest odds ratio for diabetes? Did cox.zph() suggest the proportional-hazards assumption held?

Day 5 — Reproducible Reporting & Sharing Results (Fri, Sep 25)

Inline code for a reproducible narrative

Inside a Quarto/R Markdown document (like this one!), writing a piece of inline code such as the line below directly into a sentence always reflects the current data — re-render and the numbers update themselves, with no manual copy-pasting of statistics:

The mean age was 52.1 years.
round(mean(trial$age, na.rm = TRUE), 1)
## [1] 52.1
round(sd(trial$age, na.rm = TRUE), 1)
## [1] 13.7

The mean participant age in this run was 52.1 years (SD = 13.7) — that sentence itself was generated the same way.

Exporting a Word-ready table

tbl_base <- trial |>
  group_by(treatment_arm) |>
  summarise(n = n(), mean_sbp = round(mean(sbp, na.rm = TRUE), 1), .groups = "drop")

ft <- flextable(tbl_base)
ft

treatment_arm

n

mean_sbp

Control

108

127.5

Intervention

99

129.0

save_as_docx(ft, path = "outputs/table1.docx")

Version control basics

Run in the RStudio Terminal tab, not the R console:

git init
git add .
git commit -m "Add BMI and age-group derived variables"
git remote add origin <your-github-repo-url>
git push -u origin main

CAPSTONE: build your own reproducible report

  1. Duplicate this notebook (or start a fresh .qmd/.Rmd)
  2. Swap the synthetic generate-data chunk for your own read_csv() / read_excel() import
  3. Re-run every chunk top to bottom — confirm it completes with no errors
  4. Session > Restart R, then run again — the real test of reproducibility
  5. Click Preview (or rmarkdown::render()) to produce a shareable .nb.html or .docx
cat("Notebook complete -- Day 1 through Day 5.\n")
## Notebook complete -- Day 1 through Day 5.