R for Biostatistics Exercises: 20 Real-World Practice Problems

Explanation: Table 1 in a clinical study report (CSR) summarises baseline characteristics by arm so a reviewer can confirm balance after randomization. group_by(arm) |> summarise() produces one row per arm; passing .groups = "drop" returns an ungrouped tibble so downstream code does not inherit a stray grouping. For production tables, gtsummary::tbl_summary(by = arm) adds p values and formatted output, but the raw aggregation is what regulators inspect.

Explanation: A large p value here is the desired result: it says we cannot reject independence between sex and arm, so randomization did not accidentally concentrate one sex into one arm. Yates' continuity correction defaults to on for 2 by 2 tables; for sparse cells (any expected count below 5) switch to fisher.test(tab). Imbalanced strata at baseline force the analysis to adjust covariates in the primary model, which is a regulatory complication worth avoiding.

Explanation: ICH GCP defines a protocol deviation as any departure from the approved protocol; dosing windows are a common offender. Rather than filter on absolute mg (which scales with the dose level), the percent deviation is comparable across cohorts. Always compute on abs(diff) so under and over dosing both trigger. In production, the threshold and the columns surfaced (site_id, deviation_category) come from a Statistical Analysis Plan section.

Explanation: ITT keeps every randomized subject (the gold standard for efficacy claims because it preserves randomization), while PP restricts to subjects who completed per protocol (informs adherence questions). A large arm asymmetric dropout (drug arm losing more subjects than placebo) is itself a safety signal, so this summary is the first thing a Data Safety Monitoring Board reviews. The CONSORT diagram in the publication is built from exactly this aggregation.

Explanation: Welch is the default because it does not assume equal variances, which is safer when arm sample sizes or spreads differ. The formula = y ~ group syntax is the trial reporting convention because it mirrors the model card. A 95 percent confidence interval that excludes 0 corresponds to p below 0.05 at the same alpha; report both because regulators care more about the effect size and its precision than the p value alone. For skewed endpoints, switch to wilcox.test() (see exercise 2.3).

Explanation: Pairing strips between subject variance and tests only the within subject change, which sharply boosts power for crossover and longitudinal designs. Forgetting paired = TRUE here would shrink the t statistic and inflate the confidence interval, because the unpaired test would treat the 50 observations as independent. If the differences are non normal, the paired analogue is wilcox.test(..., paired = TRUE). The reporting convention is "mean change (95 percent CI), n = X".

Explanation: Wilcoxon rank sum tests for a location shift on the ranks rather than the raw values, so heavy right tails (common for transaminases, bilirubin, troponin) do not pull the test statistic around. exact = FALSE switches to a fast normal approximation, which is the right call once either group exceeds about 20 ranks. Report a Hodges Lehmann pseudomedian (conf.int = TRUE) alongside the p value because a non parametric test still owes the reader an effect size.

Explanation: Whenever an expected cell count drops below 5, Fisher exact is preferred over chisq.test() because it conditions on the row and column marginals and computes the exact hypergeometric probability. The reported odds ratio is conditional on the margins, which is a slightly different estimand than the unconditional MLE you would get from glm(..., family = binomial); pick one and disclose it. For larger tables, switch to fisher.test(..., simulate.p.value = TRUE) to keep runtime sane.

Explanation: Surv(time, event) is a two column object recording follow up time and whether the event occurred (or the subject was censored at that time). In lung, the original coding stores 1 for censored and 2 for death, so the status == 2 predicate creates the 0/1 indicator the model needs. The reported median (310 days) is the time at which the Kaplan Meier estimate crosses 0.5; if fewer than half the subjects have died, the median is NA. For a graphical KM curve, pass the fit to plot() or survminer::ggsurvplot().

Explanation: The log rank test (Mantel Cox) is the standard test of equality between survival curves; it weights each event time equally so it has maximum power when hazards are proportional. The tiny p value here (7e-05) tells you ECOG groups differ on survival, but it does not tell you which pairs differ or by how much: for that, fit a Cox model (exercise 3.3) and read the per level hazard ratios. When proportional hazards is doubtful, switch to survdiff(..., rho = 1) for Peto Peto, which weights early events more.

Explanation: A Cox model leaves the baseline hazard unspecified and models only how covariates multiplicatively shift it, which is why hazard ratios (exp(coef)) are the reported effect size. An HR of 1.56 for ECOG means each one point worse performance score is associated with a 56 percent higher hazard of death at any moment, holding age fixed. The likelihood ratio test (preferred over Wald for small samples) tests the global null that both coefficients are zero. Report HR with 95 percent CI from confint(ex_3_3) and exp() them for the paper.

Explanation: cox.zph() regresses scaled Schoenfeld residuals on a function of time (default: KM transform) and tests the slope against zero. Per covariate p values flag whether that covariate's HR drifts; the GLOBAL row pools them. Here the global p value of 0.17 does not reject PH, so the Cox model is defensible. If PH fails, the regulator expects either a stratified Cox model (strata(ph.ecog)), a time varying coefficient (tt()), or a switch to an accelerated failure time model. A graphical companion is plot(ex_3_4).

Explanation: The random intercept (1 | Subject) gives each subject their own baseline reaction time, soaking up the obvious between subject variation (some people are just slower). The fixed slope of 10.47 ms per day of sleep restriction is the population average effect, which is the parameter a journal would report. The two variance components (Subject 37.12, Residual 30.99) say roughly half the noise is between subject and half within, which is a useful sanity check for any cluster randomized or longitudinal design.

Explanation: When comparing nested models that differ in random effect structure, refit with maximum likelihood (REML = FALSE) because REML log likelihoods are only comparable for models with identical fixed effects. The chi square of 42 on 2 df (the variance of the random slope and its covariance with the intercept) gives an absurdly small p value, telling you sleep deprivation rates vary meaningfully across people. Failure to include the random slope when it is warranted yields anti conservative standard errors on the fixed effect, the classic Type I error inflation Barr et al. flagged.

Explanation: ranef() returns each subject's BLUP (best linear unbiased predictor), which is the shrunken deviation from the population intercept. Adding it back to fixef()["(Intercept)"] reconstructs the per subject expected reaction time at Day 0. Shrinkage pulls subjects with few observations toward the population mean (a Bayesian style regularization), which is exactly why mixed models outperform per subject regressions when data is sparse. For uncertainty intervals on the BLUPs, use arm::sim() or merTools::predictInterval().

Explanation: With a single two level factor, one way ANOVA returns an F statistic that is exactly the square of the equal variance two sample t statistic; the p value here (0.06) sits at the conventional decision boundary. The borderline result is misleading because the design crosses supp with dose, and dose carries most of the variance. The honest model is a two way ANOVA that includes both factors plus their interaction (exercise 5.2), which is also the analysis pre specified in any well designed factorial experiment.

Explanation: Wrapping dose in factor() lets the model estimate a separate mean per dose rather than enforcing a linear trend; with three levels this matters because the dose effect is not perfectly linear. The supp * factor(dose) shorthand expands to main effects plus the two way interaction. The significant interaction (p = 0.022) says the supplement effect depends on dose: at 0.5 and 1.0 mg/day VC outperforms OJ; at 2.0 mg/day they converge. Always interpret interactions before main effects when both are significant.

Explanation: Logistic regression models the log odds of a binary outcome; exp(coef) converts a coefficient to an odds ratio per unit increase in the covariate. Here each additional milligram of dose multiplies the odds of seroconversion by exp(0.059) = 1.06, holding age fixed; each additional year of age divides them by exp(0.043) = 1.04. Always report odds ratios with 95 percent CI (exp(confint(ex_5_3))) and the deviance drop versus a null model as a global fit summary. For prediction calibration, use pROC::roc().

Explanation: Sample size is driven by the standardised effect size delta / sd = 0.67, the alpha, and the target power. The function returns 36.3 per arm, so the protocol rounds up to 37 per arm and adds a dropout inflation factor (e.g. divide by 1 minus expected dropout rate). For unequal arms, switch to power.t.test with n = NULL and solve over a grid, or use the pwr package which exposes pwr.t2n.test() for unbalanced designs. Always run a sensitivity sweep across plausible deltas before locking the protocol.

Explanation: Benjamini Hochberg controls the expected proportion of false discoveries among rejected hypotheses (FDR), which is a much friendlier criterion than family wise error rate when you are screening many contrasts and can tolerate a few false positives. For each rank i, BH compares p[i] against i / m * q; the largest i that satisfies the inequality defines the cut, and everything at or below it is declared significant. Bonferroni (method = "bonferroni") is strict but loses power; BH is the default for omics, GWAS, and most modern multiplicity adjustments.