# ########################################################## # Biostatistics in R — full course R scripts # https://rverseanalytics.com/courses/biostatistics # ########################################################## # ---------------------------------------------------------- # Table 1 with gtsummary # ---------------------------------------------------------- # Table 1 with gtsummary - baseline characteristics library(gtsummary) library(dplyr) # Summarise the trial cohort by treatment arm, with a p-value trial |> select(trt, age, grade, stage, marker) |> tbl_summary(by = trt, missing_text = "(missing)") |> add_p() # Add labels, an overall column, and bold row labels trial |> select(trt, age, grade, stage, marker) |> tbl_summary(by = trt, label = list(age ~ "Age (years)", marker ~ "Marker (ng/mL)")) |> add_p() |> add_overall() |> bold_labels() # ---------------------------------------------------------- # Correlation: Pearson & Spearman # ---------------------------------------------------------- # Correlation - Pearson vs Spearman # Built-in mtcars: weight vs fuel efficiency cor(mtcars$wt, mtcars$mpg) # Pearson correlation with a test and 95% CI cor.test(mtcars$wt, mtcars$mpg) # Spearman - rank based, for non-linear or non-normal data cor.test(mtcars$wt, mtcars$mpg, method = "spearman") # See the relationship plot(mtcars$wt, mtcars$mpg, pch = 19, col = "#17a2b8", xlab = "Weight", ylab = "MPG") # ---------------------------------------------------------- # Chi-square & Fisher's Test # ---------------------------------------------------------- # Chi-square & Fisher's exact test # A 2x2 table: treatment (rows) by outcome (columns) tbl <- matrix(c(90, 60, 30, 120), nrow = 2, byrow = TRUE) tbl # Chi-square test of independence chisq.test(tbl) # Expected counts under independence chisq.test(tbl)$expected # Fisher's exact test - exact, good for small samples fisher.test(tbl) # ---------------------------------------------------------- # Odds Ratio, Relative Risk & NNT # ---------------------------------------------------------- # Odds ratio, relative risk & NNT from a 2x2 table # exposed group: events / non-events exp_event <- 30; exp_none <- 70 # control group: events / non-events ctrl_event <- 15; ctrl_none <- 85 # Odds ratio and relative risk OR <- (exp_event * ctrl_none) / (exp_none * ctrl_event) RR <- (exp_event / (exp_event + exp_none)) / (ctrl_event / (ctrl_event + ctrl_none)) c(odds_ratio = OR, relative_risk = RR) # Absolute risk reduction & number needed to treat risk_exp <- exp_event / (exp_event + exp_none) risk_ctrl <- ctrl_event / (ctrl_event + ctrl_none) ARR <- risk_ctrl - risk_exp c(ARR = ARR, NNT = 1 / abs(ARR)) # ---------------------------------------------------------- # ANOVA & Post-hoc Tests # ---------------------------------------------------------- # One-way ANOVA & post-hoc # Do iris species differ in sepal length? fit <- aov(Sepal.Length ~ Species, data = iris) summary(fit) # Tukey post-hoc: which pairs differ? TukeyHSD(fit) # Visualise the groups boxplot(Sepal.Length ~ Species, data = iris, col = "#2f6fed33", border = "#1b2a4a", ylab = "Sepal length") # ---------------------------------------------------------- # Non-parametric Tests # ---------------------------------------------------------- # Non-parametric tests - no normality assumption # Mann-Whitney / Wilcoxon: two independent groups wilcox.test(mpg ~ am, data = mtcars) # Kruskal-Wallis: three or more groups kruskal.test(Sepal.Length ~ Species, data = iris) # Wilcoxon signed-rank: paired data with(sleep, wilcox.test(extra[group == 1], extra[group == 2], paired = TRUE)) # ---------------------------------------------------------- # Effect Size: Cohen's d & Hedges' g # ---------------------------------------------------------- # Effect size: Cohen's d & Hedges' g library(effectsize) # How big is the mpg gap between gearboxes? cohens_d(mpg ~ am, data = mtcars) # Hedges' g - bias-corrected for small samples hedges_g(mpg ~ am, data = mtcars) # Interpret the magnitude interpret_cohens_d(1.48) # ---------------------------------------------------------- # Multiple Linear Regression # ---------------------------------------------------------- # Multiple linear regression library(gtsummary) # Predict mpg from weight, horsepower and gearbox fit <- lm(mpg ~ wt + hp + am, data = mtcars) summary(fit) # A tidy regression table tbl_regression(fit) # ---------------------------------------------------------- # Logistic Regression & Odds Ratios # ---------------------------------------------------------- # Logistic regression - odds ratios & ROC library(gtsummary) library(pROC) # Fit a logistic model for tumour response d <- na.omit(trial[, c("response", "age", "grade", "stage")]) fit <- glm(response ~ age + grade + stage, data = d, family = binomial) # Adjusted odds ratios with 95% CI tbl_regression(fit, exponentiate = TRUE) |> bold_p() # ROC curve and AUC - how well the model discriminates roc_obj <- roc(d$response, fitted(fit), quiet = TRUE) auc(roc_obj) plot(roc_obj, col = "#2f6fed", lwd = 3, legacy.axes = TRUE) # ---------------------------------------------------------- # Kaplan-Meier Survival Curves # ---------------------------------------------------------- # Kaplan-Meier survival curves library(survival) # Survival by sex (lung cancer data) fit <- survfit(Surv(time, status) ~ sex, data = lung) fit # The survival curve plot(fit, col = c("#2f6fed", "#17a2b8"), lwd = 3, xlab = "Days", ylab = "Survival probability") legend("topright", c("Male", "Female"), col = c("#2f6fed", "#17a2b8"), lwd = 3) # Log-rank test: do the groups differ? survdiff(Surv(time, status) ~ sex, data = lung) # ---------------------------------------------------------- # Cox Proportional Hazards # ---------------------------------------------------------- # Cox proportional hazards library(survival) library(gtsummary) # Fit a Cox model (lung cancer data) cx <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung) # Hazard ratios with 95% CI tbl_regression(cx, exponentiate = TRUE) # Check the proportional-hazards assumption cox.zph(cx) # ---------------------------------------------------------- # Diagnostic Test Accuracy # ---------------------------------------------------------- # Diagnostic test accuracy from a 2x2 table # test positive: TP true, FP false ; test negative: FN, TN TP <- 90; FP <- 10 FN <- 20; TN <- 80 # Sensitivity & specificity sensitivity <- TP / (TP + FN) specificity <- TN / (TN + FP) c(sensitivity = sensitivity, specificity = specificity) # Predictive values PPV <- TP / (TP + FP) NPV <- TN / (TN + FN) c(PPV = PPV, NPV = NPV) # Likelihood ratios c(LR_positive = sensitivity / (1 - specificity), LR_negative = (1 - sensitivity) / specificity) # ---------------------------------------------------------- # Cohen's Kappa: Rater Agreement # ---------------------------------------------------------- # Inter-rater agreement: Cohen's kappa (base R) # Two raters classify 50 items; their agreement table: tab <- matrix(c(20, 5, 3, 22), nrow = 2, byrow = TRUE, dimnames = list(Rater1 = c("+", "-"), Rater2 = c("+", "-"))) tab # Observed vs chance-expected agreement po <- sum(diag(tab)) / sum(tab) pe <- sum(rowSums(tab) * colSums(tab)) / sum(tab)^2 c(observed = po, expected = pe) # Cohen's kappa - agreement beyond chance kappa <- (po - pe) / (1 - pe) kappa # ---------------------------------------------------------- # Cronbach's Alpha: Reliability # ---------------------------------------------------------- # Cronbach's alpha - scale reliability library(psych) # Five agreeableness items from the bfi dataset items <- bfi[, c("A1", "A2", "A3", "A4", "A5")] # Cronbach's alpha (A1 is reverse-keyed) a <- alpha(items, check.keys = TRUE) a$total # ---------------------------------------------------------- # Sample Size & Power # ---------------------------------------------------------- # Sample size & power for a two-group t-test # n per group for 80% power (delta = 5, sd = 10) power.t.test(delta = 5, sd = 10, power = 0.80, sig.level = 0.05) # Power for a fixed n = 30 per group power.t.test(n = 30, delta = 5, sd = 10, sig.level = 0.05) # With the pwr package: effect size d = 0.5 library(pwr) pwr.t.test(d = 0.5, power = 0.80, sig.level = 0.05) # Two proportions: 40% vs 60% pwr.2p.test(h = ES.h(0.6, 0.4), power = 0.80)