V The V Lab

Design of Experiments: Principles, Analysis, and Applications

From randomization, replication, and blocking through factorial effects, randomization inference, power, diagnostics, and responsible reporting

The V Lab Home Biostatistics Foundations Medical Statistical Tests Linear Regression Logistic Regression Design of Experiments Clinical Trials Poisson and Zero-Inflated Models Multilevel Modelling Longitudinal Data Analysis Causal Inference Propensity Score Matching Survival Analysis AFT and Cox PH Meta-Analysis 中文版

AudienceLearners in public health, epidemiology, medicine, and health data science

Study timeApproximately 180–240 minutes

PrerequisitesMeans, variance, confidence intervals, linear models, and basic R

About the simulated experiment This module follows a fixed-seed, simulated 2 × 2 factorial trial in 12 community clinics. Within every clinic, four participants receive each combination of health coaching and home blood-pressure monitoring, for 192 participants total. The outcome is 12-week reduction in systolic blood pressure (mmHg), where a larger value is better. No real people or clinical records are represented.

How to use this tutorial

Define the causal question and experimental unit → choose factors, levels, controls, and blocks → randomize reproducibly → preserve the assignment structure in analysis → estimate factorial contrasts → diagnose assumptions and attrition → report estimands, uncertainty, and design limitations.

Learning objectives

After completing this tutorial, you should be able to:

• identify the experimental unit, observational unit, factor, level, treatment combination, response, block, replicate, and estimand;
• explain why randomization, replication, local control, allocation concealment, blinding, and preregistration solve different problems;
• distinguish a completely randomized design (CRD) from a randomized complete block design (RCBD);
• design and audit a balanced 2 × 2 factorial experiment;
• distinguish simple effects, marginal main effects, and interaction effects;
• fit and interpret lm(bp_reduction ~ clinic + coaching * home_monitoring);
• obtain cell means and valid linear-combination confidence intervals averaged over all clinics;
• use prespecified contrasts, multiplicity-aware post hoc comparisons, residual diagnostics, and restricted randomization inference;
• plan sample size with an explicit effect, variance, analysis, and simulation model;
• preserve intention-to-treat logic under nonadherence and address missing outcomes transparently; and
• recognize when cluster, split-plot, crossover, Latin-square, repeated-measures, binary, or count designs need a different analysis.

1 Start with the scientific question

1.1 Treatments, units, outcomes, and estimands

An experiment deliberately assigns one or more interventions and observes their consequences. Design begins before a statistical test is selected.

Term	Meaning in this trial	Question to ask
Experimental unit	Participant receiving an assigned combination	What entity was independently randomized?
Observational unit	Participant whose 12-week blood-pressure change is measured	At what level is the response recorded?
Factor	Coaching; home monitoring	What intervention dimension is manipulated?
Level	No or Yes for each factor	Which versions are compared?
Treatment combination	Control, coaching only, monitoring only, or both	Which joint factor levels can be assigned?
Block	Clinic	Which known nuisance source is controlled locally?
Replicate	A distinct randomized participant at a combination	Is there independent information about experimental error?
Response	Systolic BP reduction at 12 weeks, in mmHg	Is direction, timing, and measurement protocol prespecified?
Estimand	For example, marginal effect of coaching averaged over monitoring levels and clinics	Exactly which contrast in potential outcomes is targeted?

The experimental unit follows randomization, not the row count. Here, participants are randomized within clinics, so 192 participants are experimental units and clinics are blocks. If a whole clinic received one policy, the clinic—not each patient—would be the experimental unit.

Let $Y_i(a,b)$ be participant $i$’s potential BP reduction under coaching level $a\in\{0,1\}$ and monitoring level $b\in\{0,1\}$. One finite-sample marginal coaching estimand is

$$\tau_A=\frac{1}{2}\left\{E[Y(1,0)-Y(0,0)]+E[Y(1,1)-Y(0,1)]\right\}.$$

The difference-in-differences interaction estimand is

$$\tau_{AB}=E[Y(1,1)-Y(0,1)]-E[Y(1,0)-Y(0,0)].$$

These definitions specify how factor levels are averaged. They are clearer than saying only “the treatment effect.” Generalization from these trial participants to another population is a separate question from internal causal identification.

A nonsignificant result is not evidence of no effect A large $p$-value can coexist with effects that are clinically important in either direction. Report the prespecified effect estimate, confidence interval, outcome scale, and compatibility with clinically meaningful values. If equivalence or noninferiority is the goal, design and analyze that question explicitly.

1.2 Six design principles with distinct jobs

1. Randomization makes assignment probabilities known and protects against systematic baseline differences in expectation. It supports design-based inference when the actual randomization is respected.
2. Replication supplies independent experimental information and permits estimation of random variation. Repeated measurements on one unit are subsamples, not independent treatment replicates.
3. Blocking or local control compares treatments within relatively homogeneous groups, often improving precision and guaranteeing local balance.
4. Control conditions define the counterfactual comparison: usual care, placebo, attention control, active comparator, or another policy. The choice determines the estimand.
5. Allocation concealment and blinding address different biases. Concealment prevents foreknowledge before assignment; blinding limits postassignment differences in care, behavior, measurement, or assessment when feasible.
6. Preregistration and an analysis plan distinguish confirmatory outcomes and contrasts from exploratory work, document stopping rules and exclusions, and reduce undisclosed analytical flexibility.

Randomization does not guarantee perfect baseline balance in one realized trial, cure attrition, ensure adherence, prevent measurement bias, or make an unrepresentative sample representative. Those problems need their own design and analysis protections.

2 The 2 × 2 clinic-blocked experiment

2.1 Data-generating mechanism and audit

Each clinic contains all four combinations, with four independently randomized participants per combination. Before looking at outcomes, audit unique identifiers, factor levels, allocation counts, and missing values.

design_audit <- data.frame(
  Item = c(
    "Participants", "Clinics", "Participants per clinic",
    "Participants per clinic-combination", "Missing outcomes",
    "Duplicate clinic-participant identifiers"
  ),
  Value = c(
    nrow(doe_data), nlevels(doe_data[["clinic"]]),
    paste(range(as.vector(table(doe_data[["clinic"]]))), collapse = " to "),
    paste(range(as.vector(xtabs(~ clinic + treatment_combination, doe_data))),
          collapse = " to "),
    sum(is.na(doe_data[["bp_reduction"]])),
    sum(duplicated(doe_data[c("clinic", "participant_in_clinic")]))
  )
)
knitr::kable(design_audit, caption = "Pre-outcome structural audit of the simulated experiment")

Pre-outcome structural audit of the simulated experiment

Item	Value
Participants	192
Clinics	12
Participants per clinic	16 to 16
Participants per clinic-combination	4 to 4
Missing outcomes	0
Duplicate clinic-participant identifiers	0

The simulation has true clinic effects with standard deviation 2.2 mmHg and individual errors with standard deviation 4.5 mmHg. Its four population cell means, before clinic shifts, are 4.0, 7.0, 6.2, and 12.0 mmHg. Therefore:

• the coaching simple effect is 3.0 mmHg without monitoring and 5.8 mmHg with monitoring;
• the monitoring simple effect is 2.2 mmHg without coaching and 5.0 mmHg with coaching;
• the marginal coaching effect is $(3.0+5.8)/2=4.4$ mmHg;
• the marginal monitoring effect is $(2.2+5.0)/2=3.6$ mmHg; and
• the difference-in-differences interaction is $5.8-3.0=2.8$ mmHg.

These are simulation truths, not estimates available in a real study.

2.2 A reproducible randomization schedule

A randomization list is part of the design record. It should be generated by an authorized reproducible procedure, secured before enrollment, and paired with identifiers only through the allocation system. Never reconstruct a sequence after seeing outcomes.

schedule_preview <- doe_data[
  doe_data[["clinic"]] %in% clinic_levels[1:2],
  c("clinic", "participant_in_clinic", "treatment_combination")
]
knitr::kable(schedule_preview,
  caption = "First two clinic-specific assignment schedules (illustrative only)")

First two clinic-specific assignment schedules (illustrative only)

clinic	participant_in_clinic	treatment_combination
Clinic 01	1	Both
Clinic 01	2	Monitoring only
Clinic 01	3	Coaching only
Clinic 01	4	Monitoring only
Clinic 01	5	Control
Clinic 01	6	Monitoring only
Clinic 01	7	Coaching only
Clinic 01	8	Both
Clinic 01	9	Coaching only
Clinic 01	10	Coaching only
Clinic 01	11	Control
Clinic 01	12	Control
Clinic 01	13	Both
Clinic 01	14	Both
Clinic 01	15	Control
Clinic 01	16	Monitoring only
Clinic 02	1	Coaching only
Clinic 02	2	Coaching only
Clinic 02	3	Coaching only
Clinic 02	4	Both
Clinic 02	5	Monitoring only
Clinic 02	6	Coaching only
Clinic 02	7	Both
Clinic 02	8	Control
Clinic 02	9	Both
Clinic 02	10	Monitoring only
Clinic 02	11	Control
Clinic 02	12	Control
Clinic 02	13	Control
Clinic 02	14	Both
Clinic 02	15	Monitoring only
Clinic 02	16	Monitoring only

allocation_matrix <- matrix(doe_data[["combination"]], nrow = n_clinics, byrow = TRUE)
old_par <- par(mar = c(5, 6, 4, 1))
image(
  x = seq_len(16), y = seq_len(n_clinics), z = t(allocation_matrix),
  col = c(palette_doe["gray"], palette_doe["blue"],
          palette_doe["orange"], palette_doe["teal"]),
  breaks = seq(0.5, 4.5, by = 1), axes = FALSE,
  xlab = "Enrollment position within clinic", ylab = "Clinic",
  main = "Randomized complete block allocation"
)
axis(1, at = c(1, 4, 8, 12, 16))
axis(2, at = seq_len(n_clinics), labels = clinic_levels, las = 1, cex.axis = 0.7)
legend(
  "bottom", inset = 0.01, horiz = TRUE, bty = "n",
  legend = levels(doe_data[["treatment_combination"]]),
  fill = c(palette_doe["gray"], palette_doe["blue"],
           palette_doe["orange"], palette_doe["teal"]), cex = 0.8
)

Clinic-specific randomized treatment schedules. Every row is a clinic, every column is an enrollment position, and color indicates one of the four treatment combinations. The order varies while each clinic retains four assignments to every combination.

par(old_par)

3 CRD, blocking, and factorial structure

3.1 Completely randomized design

In a completely randomized design, treatment combinations are assigned across all eligible experimental units without a blocking restriction. With 192 participants and four arms, a balanced CRD could be generated as follows.

set.seed(314159)
crd_assignment <- sample(rep(
  c("Control", "Coaching only", "Monitoring only", "Both"), each = 48
))
knitr::kable(as.data.frame(table(crd_assignment)),
  col.names = c("Treatment combination", "Assigned participants"),
  caption = "Balanced completely randomized allocation example")

Balanced completely randomized allocation example

Treatment combination	Assigned participants
Both	48
Coaching only	48
Control	48
Monitoring only	48

A CRD is simple and flexible. It can be efficient when units are homogeneous or no strong preassignment prognostic factor is available. In a multiclinic trial, however, chance might place disproportionate combinations in clinics with systematically different populations or measurement practices.

3.2 Randomized complete block design

In an RCBD, each block contains every treatment combination. Here, each clinic independently randomizes four participants to each combination. The analysis compares combinations after accounting for clinic, and the randomization test must permute only within clinics.

allocation_by_clinic <- xtabs(~ clinic + treatment_combination, data = doe_data)
knitr::kable(allocation_by_clinic,
  caption = "Treatment-combination counts within all 12 complete clinic blocks")

Treatment-combination counts within all 12 complete clinic blocks

	Control	Coaching only	Monitoring only	Both
Clinic 01	4	4	4	4
Clinic 02	4	4	4	4
Clinic 03	4	4	4	4
Clinic 04	4	4	4	4
Clinic 05	4	4	4	4
Clinic 06	4	4	4	4
Clinic 07	4	4	4	4
Clinic 08	4	4	4	4
Clinic 09	4	4	4	4
Clinic 10	4	4	4	4
Clinic 11	4	4	4	4
Clinic 12	4	4	4	4

Blocking is most useful when blocks explain outcome variation and every planned treatment comparison occurs within blocks. Blocks should be defined by preassignment information. Excessively many tiny strata can make implementation fragile; dynamic allocation or covariate-adaptive procedures require their own documented inferential methods.

3.3 Why a factorial experiment?

A 2 × 2 factorial design learns about two factors and their interaction in the same experiment. With equal allocation, each coaching main-effect comparison uses all 192 participants: coaching is present in half and absent in half, averaged across monitoring. This can be much more efficient than running two unrelated trials—provided the joint interventions are feasible and the interaction is scientifically meaningful.

For cell means $\mu_{00},\mu_{10},\mu_{01},\mu_{11}$ in the order control, coaching only, monitoring only, and both:

$$\begin{aligned} \text{Marginal coaching effect} &= \tfrac12(-\mu_{00}+\mu_{10}-\mu_{01}+\mu_{11}),\\ \text{Marginal monitoring effect} &= \tfrac12(-\mu_{00}-\mu_{10}+\mu_{01}+\mu_{11}),\\ \text{Interaction} &= \mu_{00}-\mu_{10}-\mu_{01}+\mu_{11}. \end{aligned}$$

An interaction means the effect of one factor differs across levels of the other. It does not automatically imply a biological mechanism, and its scale matters: additivity on the mmHg scale is a different assumption from additivity on a log-risk or log-rate scale.

Coding determines coefficient labels, not the estimand With R’s default 0/1 treatment coding, the coachingYes coefficient is coaching’s simple effect when monitoring is “No,” and home_monitoringYes is monitoring’s simple effect when coaching is “No.” The interaction coefficient is the difference in differences. With $-1/+1$ effect coding in $\mu=\gamma_0+\gamma_A A+\gamma_B B+\gamma_{AB}AB$, the marginal effects are $2\gamma_A$ and $2\gamma_B$, while the interaction difference in differences is $4\gamma_{AB}$. Always state the contrast itself.

4 Explore outcomes without breaking the design

4.1 Cell summaries

Descriptive statistics should preserve the randomized combinations and blocks. They are not substitutes for uncertainty estimates, and observed baseline or outcome imbalances should not trigger undocumented changes to a prespecified analysis.

cell_n <- aggregate(bp_reduction ~ coaching + home_monitoring, doe_data, length)
cell_mean <- aggregate(bp_reduction ~ coaching + home_monitoring, doe_data, mean)
cell_sd <- aggregate(bp_reduction ~ coaching + home_monitoring, doe_data, sd)
names(cell_n)[3] <- "n"
names(cell_mean)[3] <- "Mean"
names(cell_sd)[3] <- "SD"
cell_summary <- Reduce(
  function(x, y) merge(x, y, by = c("coaching", "home_monitoring")),
  list(cell_n, cell_mean, cell_sd)
)
knitr::kable(cell_summary, digits = 2,
  caption = "Observed BP-reduction summaries by randomized factorial cell")

Observed BP-reduction summaries by randomized factorial cell

coaching	home_monitoring	n	Mean	SD
No	No	48	2.96	4.72
No	Yes	48	6.18	5.55
Yes	No	48	5.77	6.23
Yes	Yes	48	11.24	5.08

boxplot(
  bp_reduction ~ treatment_combination, data = doe_data,
  col = c("#E8ECEF", "#C9E1F2", "#F8E2B5", "#B9DFDB"), border = palette_doe["navy"],
  ylab = "12-week systolic BP reduction (mmHg)", xlab = "",
  main = "Outcome distributions by treatment combination", las = 1
)
stripchart(
  bp_reduction ~ treatment_combination, data = doe_data,
  vertical = TRUE, method = "jitter", add = TRUE, pch = 16,
  col = grDevices::adjustcolor(palette_doe["navy"], alpha.f = 0.35)
)
abline(h = 0, lty = 3, col = palette_doe["gray"])

Observed 12-week systolic blood-pressure reduction by randomized treatment combination. Each box summarizes 48 participants across all 12 clinics; points show individual outcomes and should not be mistaken for independent clinic-level assignments.

4.2 Interaction plot

Parallel traces suggest additivity on the plotted scale; nonparallel traces suggest interaction. Sampling noise can also create nonparallel lines, so the plot should accompany a contrast estimate and interval.

with(doe_data, interaction.plot(
  x.factor = home_monitoring, trace.factor = coaching,
  response = bp_reduction, fun = mean, type = "b", pch = c(1, 19), lwd = 2,
  col = c(palette_doe["orange"], palette_doe["teal"]),
  xlab = "Home monitoring", ylab = "Mean BP reduction (mmHg)",
  trace.label = "Coaching", main = "Observed factorial pattern"
))

Mean BP reduction for coaching and no coaching across home-monitoring levels. Nonparallel lines visualize the positive interaction: the coaching difference is larger when home monitoring is also provided.

5 Analyze the randomized complete block factorial design

5.1 Fit the prespecified linear model

For this balanced design with a continuous outcome, use fixed clinic-block indicators plus both factor main terms and their interaction:

$$Y_{ijk}=\beta_0+\alpha_j+\beta_A A_{ijk}+\beta_B B_{ijk} +\beta_{AB}A_{ijk}B_{ijk}+\varepsilon_{ijk}.$$

The clinic coefficients remove between-clinic level shifts. The model assumes independent experimental-unit errors with a common variance and an adequate additive clinic structure. Randomization supplies the causal comparison; the outcome model supplies a convenient estimator and standard error under its assumptions.

rcbd_fit <- lm(
  bp_reduction ~ clinic + coaching * home_monitoring,
  data = doe_data
)
treatment_terms <- c(
  "coachingYes", "home_monitoringYes",
  "coachingYes:home_monitoringYes"
)
coefficient_table <- coef(summary(rcbd_fit))[treatment_terms, , drop = FALSE]
coefficient_table <- data.frame(
  Term = c(
    "Coaching simple effect when monitoring = No",
    "Monitoring simple effect when coaching = No",
    "Interaction: difference in differences"
  ),
  Estimate = coefficient_table[, "Estimate"],
  SE = coefficient_table[, "Std. Error"],
  t = coefficient_table[, "t value"],
  p_value = format_p(coefficient_table[, "Pr(>|t|)"]), row.names = NULL
)
knitr::kable(coefficient_table, digits = 3,
  caption = "Treatment-coded coefficients from the clinic-blocked factorial model")

Treatment-coded coefficients from the clinic-blocked factorial model

Term	Estimate	SE	t	p_value
Coaching simple effect when monitoring = No	2.81	0.96	2.92	0.00392
Monitoring simple effect when coaching = No	3.21	0.96	3.35	< 0.001
Interaction: difference in differences	2.25	1.36	1.66	0.09860

The first two coefficient rows are simple effects at reference levels, not marginal main effects. Their scientific meaning changes if the reference level changes. The next section estimates named contrasts that do not depend on a reader remembering the coding convention.

5.2 Balanced-design ANOVA decomposition

aov_fit <- aov(
  bp_reduction ~ clinic + coaching * home_monitoring,
  data = doe_data
)
anova_table <- anova(aov_fit)
anova_display <- data.frame(
  Term = rownames(anova_table), anova_table,
  row.names = NULL, check.names = FALSE
)
anova_display[["Pr(>F)"]] <- format_p(anova_display[["Pr(>F)"]])
knitr::kable(anova_display, digits = 3,
  caption = "ANOVA decomposition for clinic blocks and the 2 × 2 factorial treatment")

ANOVA decomposition for clinic blocks and the 2 × 2 factorial treatment

Term	Df	Sum Sq	Mean Sq	F value	Pr(>F)
clinic	11	1617	147.0	6.64	<0.001
coaching	1	743	742.7	33.57	<0.001
home_monitoring	1	904	904.5	40.88	<0.001
coaching:home_monitoring	1	61	61.0	2.76	0.0986
Residuals	177	3917	22.1	NA	NA

Because every clinic contains equal replication of all four combinations, the clinic, marginal factor, and interaction comparisons are orthogonal in this complete balanced design. Consequently, the treatment decomposition has a clean design-based interpretation. In an unbalanced dataset—especially after missing outcomes—sequential sums of squares can depend on model order. Return to estimand-based cell-mean contrasts rather than selecting a “Type” of sums of squares mechanically.

5.3 Estimate cell means averaged over all clinics

An intercept describes the reference clinic and control cell, which is usually not the desired population summary. Construct one design-matrix row for every clinic at each treatment combination, average those 12 rows, and propagate the full covariance matrix. This standardizes every cell to the same empirical distribution of clinics.

cell_profiles <- data.frame(
  coaching = factor(c("No", "Yes", "No", "Yes"), levels = c("No", "Yes")),
  home_monitoring = factor(c("No", "No", "Yes", "Yes"), levels = c("No", "Yes")),
  Cell = c("Control", "Coaching only", "Monitoring only", "Both")
)

average_model_row <- function(coaching_value, monitoring_value) {
  profile <- data.frame(
    clinic = factor(clinic_levels, levels = clinic_levels),
    coaching = factor(rep(coaching_value, n_clinics), levels = c("No", "Yes")),
    home_monitoring = factor(rep(monitoring_value, n_clinics), levels = c("No", "Yes"))
  )
  colMeans(model.matrix(delete.response(terms(rcbd_fit)), data = profile))
}

cell_X <- do.call(rbind, lapply(seq_len(nrow(cell_profiles)), function(j) {
  average_model_row(
    as.character(cell_profiles[["coaching"]][j]),
    as.character(cell_profiles[["home_monitoring"]][j])
  )
}))

cell_estimate <- drop(cell_X %*% coef(rcbd_fit))
cell_se <- sqrt(diag(cell_X %*% vcov(rcbd_fit) %*% t(cell_X)))
t_critical <- qt(0.975, df = df.residual(rcbd_fit))
cell_results <- data.frame(
  Cell = cell_profiles[["Cell"]],
  Estimated_mean = cell_estimate,
  SE = cell_se,
  CI_low = cell_estimate - t_critical * cell_se,
  CI_high = cell_estimate + t_critical * cell_se
)
knitr::kable(cell_results, digits = 2,
  caption = "Model-based cell means standardized equally across all 12 clinics")

Model-based cell means standardized equally across all 12 clinics

Cell	Estimated_mean	SE	CI_low	CI_high
Control	2.96	0.68	1.62	4.30
Coaching only	5.77	0.68	4.43	7.11
Monitoring only	6.18	0.68	4.84	7.52
Both	11.24	0.68	9.90	12.58

The intervals describe uncertainty in mean outcomes under this model. They are not prediction intervals for an individual, which must include individual residual variation, and they do not quantify transport uncertainty to clinics outside the study.

5.4 Prespecified marginal and interaction contrasts

Any scientifically meaningful factorial estimand can be written $L\mu$, where $\mu=(\mu_{00},\mu_{10},\mu_{01},\mu_{11})^T$. We use the model-matrix rows above so estimates, standard errors, and covariances all remain aligned.

contrast_weights <- rbind(
  "Marginal coaching effect" = c(-0.5, 0.5, -0.5, 0.5),
  "Marginal monitoring effect" = c(-0.5, -0.5, 0.5, 0.5),
  "Interaction (difference in differences)" = c(1, -1, -1, 1),
  "Coaching effect when monitoring = No" = c(-1, 1, 0, 0),
  "Coaching effect when monitoring = Yes" = c(0, 0, -1, 1)
)
contrast_X <- contrast_weights %*% cell_X
contrast_estimate <- drop(contrast_X %*% coef(rcbd_fit))
contrast_se <- sqrt(diag(contrast_X %*% vcov(rcbd_fit) %*% t(contrast_X)))
contrast_t <- contrast_estimate / contrast_se
contrast_results <- data.frame(
  Contrast = rownames(contrast_weights),
  Estimate = contrast_estimate,
  SE = contrast_se,
  CI_low = contrast_estimate - t_critical * contrast_se,
  CI_high = contrast_estimate + t_critical * contrast_se,
  nominal_unadjusted_p = format_p(
    2 * pt(abs(contrast_t), df = df.residual(rcbd_fit), lower.tail = FALSE)
  ),
  row.names = NULL
)
knitr::kable(
  contrast_results, digits = 3,
  caption = "Prespecified factorial effects and simple effects with 95% confidence intervals; p-values are nominal and must follow the prespecified multiplicity plan"
)

Prespecified factorial effects and simple effects with 95% confidence intervals; p-values are nominal and must follow the prespecified multiplicity plan

Contrast	Estimate	SE	CI_low	CI_high	nominal_unadjusted_p
Marginal coaching effect	3.93	0.679	2.594	5.27	< 0.001
Marginal monitoring effect	4.34	0.679	3.001	5.68	< 0.001
Interaction (difference in differences)	2.25	1.358	-0.425	4.93	0.09860
Coaching effect when monitoring = No	2.81	0.960	0.911	4.70	0.00392
Coaching effect when monitoring = Yes	5.06	0.960	3.166	6.96	< 0.001

The marginal coaching estimate averages coaching’s effect equally across monitoring levels, as prespecified. If the target population would use monitoring at another prevalence, replace 0.5/0.5 with justified target weights and label the estimand accordingly. When interaction is substantial, lead with cell means and simple effects; a single marginal effect can hide decision-relevant heterogeneity.

Interpretation in the simulated experiment The estimated marginal coaching effect is 3.93 mmHg (95% CI 2.59 to 5.27). The estimated interaction is 2.25 mmHg (95% CI -0.43 to 4.93). These estimates quantify assigned-intervention contrasts under the simulated design. Their intervals should be read as ranges of values compatible with the model and data, not as binary declarations that an effect exists or does not exist.

5.5 Multiple comparisons and Tukey intervals

Prespecified primary contrasts should drive sample size and interpretation. If all six pairwise cell comparisons are explored, simultaneous coverage matters. TukeyHSD() provides familywise-error-controlled intervals for pairwise comparisons in this balanced Gaussian ANOVA; it does not decide which comparisons are scientifically important.

tukey_cells <- TukeyHSD(aov_fit, which = "coaching:home_monitoring")[[1]]
tukey_display <- data.frame(
  Comparison = rownames(tukey_cells), tukey_cells,
  row.names = NULL, check.names = FALSE
)
tukey_display[["p adj"]] <- format_p(tukey_display[["p adj"]])
knitr::kable(tukey_display, digits = 3,
  caption = "Tukey-adjusted pairwise comparisons among the four factorial cells")

Tukey-adjusted pairwise comparisons among the four factorial cells

Comparison	diff	lwr	upr	p adj
Yes:No-No:No	2.806	0.316	5.30	0.02031
No:Yes-No:No	3.214	0.723	5.70	0.00547
Yes:Yes-No:No	8.275	5.784	10.77	< 0.001
No:Yes-Yes:No	0.407	-2.083	2.90	0.97430
Yes:Yes-Yes:No	5.468	2.978	7.96	< 0.001
Yes:Yes-No:Yes	5.061	2.571	7.55	< 0.001

Do not use a global ANOVA as a gatekeeper that forbids reporting a prespecified contrast. Conversely, do not call one of many data-selected comparisons “confirmatory.” Report the family of hypotheses, adjustment, and whether an analysis was planned or exploratory.

6 What blocking buys

6.1 Compare blocked and unblocked analyses

Ignoring clinic does not undo randomization, but it leaves systematic clinic variation in the residual and can reduce precision. Because treatment is balanced within every clinic, adjustment is prespecified and cannot be driven by observed outcome imbalance.

crd_fit <- lm(bp_reduction ~ coaching * home_monitoring, data = doe_data)

linear_combo_se <- function(model, weights) {
  coefficient_weights <- setNames(rep(0, length(coef(model))), names(coef(model)))
  coefficient_weights[names(weights)] <- weights
  sqrt(drop(t(coefficient_weights) %*% vcov(model) %*% coefficient_weights))
}

marginal_coaching_weights <- c(
  coachingYes = 1,
  "coachingYes:home_monitoringYes" = 0.5
)
se_coaching_crd <- linear_combo_se(crd_fit, marginal_coaching_weights)
se_coaching_rcbd <- linear_combo_se(rcbd_fit, marginal_coaching_weights)

blocking_comparison <- data.frame(
  Analysis = c("Ignore clinic blocks", "Include clinic blocks"),
  Residual_MSE = c(deviance(crd_fit) / df.residual(crd_fit),
                   deviance(rcbd_fit) / df.residual(rcbd_fit)),
  Marginal_coaching_SE = c(se_coaching_crd, se_coaching_rcbd)
)
knitr::kable(blocking_comparison, digits = 3,
  caption = "Residual variation and coaching-effect precision with and without clinic blocks")

Residual variation and coaching-effect precision with and without clinic blocks

Analysis	Residual_MSE	Marginal_coaching_SE
Ignore clinic blocks	29.4	0.783
Include clinic blocks	22.1	0.679

The estimated relative efficiency from residual mean squares is 1.33: under these assumptions, the unblocked analysis would need roughly that multiplier of information to attain comparable residual precision. This is a descriptive comparison for the simulated design, not a universal property of blocking.

clinic_means <- aggregate(bp_reduction ~ clinic, data = doe_data, FUN = mean)
dotchart(
  clinic_means[["bp_reduction"]], labels = clinic_means[["clinic"]],
  pch = 19, color = palette_doe["teal"],
  xlab = "Clinic mean BP reduction (mmHg)",
  main = "Outcome differences across clinic blocks"
)
abline(v = mean(doe_data[["bp_reduction"]]), lty = 2, lwd = 2,
       col = palette_doe["orange"])

Mean BP reduction by clinic, with an overall mean reference line. Between-clinic variation is a nuisance source absorbed by clinic blocks, leaving treatment comparisons to be made locally within every clinic.

7 Model diagnostics and robust reasoning

7.1 Residual checks

The linear-model standard errors assume independent, approximately homoscedastic errors and a correctly specified mean. Randomization protects the assignment mechanism, but it does not make every model-based standard error valid. Inspect residual patterns, tails, influence, and the level at which dependence can arise.

standardized_residuals <- rstandard(rcbd_fit)
old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 2.5, 1))
plot(
  fitted(rcbd_fit), resid(rcbd_fit), pch = 19,
  col = grDevices::adjustcolor(palette_doe["teal"], alpha.f = 0.55),
  xlab = "Fitted value", ylab = "Residual", main = "Residuals vs fitted"
)
abline(h = 0, lty = 2, col = palette_doe["gray"])
qqnorm(standardized_residuals, pch = 19, col = palette_doe["blue"],
       main = "Normal Q-Q")
qqline(standardized_residuals, col = palette_doe["orange"], lwd = 2)
plot(
  fitted(rcbd_fit), sqrt(abs(standardized_residuals)), pch = 19,
  col = grDevices::adjustcolor(palette_doe["vermillion"], alpha.f = 0.55),
  xlab = "Fitted value", ylab = expression(sqrt("|standardized residual|")),
  main = "Scale-location"
)
plot(
  cooks.distance(rcbd_fit), type = "h", col = palette_doe["navy"],
  xlab = "Participant row", ylab = "Cook's distance", main = "Influence"
)

Four diagnostics for the clinic-blocked factorial linear model: residuals versus fitted values, a normal quantile plot, scale-location behavior, and Cook’s distance. They assess mean structure, variance, tail behavior, and influential experimental units.

par(old_par)

With moderate balanced cells, mean contrasts can be reasonably robust to mild nonnormality, but heavy tails, unequal variances, protocol errors, or informative missingness deserve sensitivity analyses. A transformation changes the effect scale. Heteroskedasticity-consistent or cluster-robust standard errors can be useful in other settings, but they do not repair the wrong experimental-unit definition, and reliable cluster-robust inference needs enough independent clusters.

7.2 Restricted randomization inference

A randomization test can analyze the treatment assignment mechanism directly. An exact test of the global sharp null would enumerate every allowed assignment while keeping outcomes and clinics fixed and jointly permuting the four-level treatment-combination labels within each clinic. Below, we approximate that distribution with 999 Monte Carlo draws and a plus-one $p$-value. We prespecify two statistics: the interaction difference in differences and a three-degree-of-freedom omnibus treatment $F$ statistic.

omnibus_treatment_f <- function(data) {
  reduced_fit <- lm(bp_reduction ~ clinic, data = data)
  full_fit <- lm(bp_reduction ~ clinic + treatment_combination, data = data)
  unname(anova(reduced_fit, full_fit)[["F"]][2])
}

interaction_statistic <- function(labels, outcome) {
  means <- tapply(outcome, factor(labels, levels = 1:4), mean)
  unname((means[4] - means[3]) - (means[2] - means[1]))
}

observed_f <- omnibus_treatment_f(doe_data)
observed_interaction <- interaction_statistic(
  doe_data[["combination"]], doe_data[["bp_reduction"]]
)
set.seed(20260915)
B_permutations <- 999
clinic_indices <- split(seq_len(nrow(doe_data)), doe_data[["clinic"]])

permuted_f <- replicate(B_permutations, {
  permuted_data <- doe_data
  permuted_labels <- as.character(doe_data[["treatment_combination"]])
  for (indices in clinic_indices) {
    permuted_labels[indices] <- sample(permuted_labels[indices], replace = FALSE)
  }
  permuted_data[["treatment_combination"]] <- factor(
    permuted_labels, levels = levels(doe_data[["treatment_combination"]])
  )
  omnibus_treatment_f(permuted_data)
})

# Reset the seed so the focused statistic uses exactly the same allowed
# assignment sequence as the omnibus statistic and the Chinese companion page.
set.seed(20260915)
permuted_interaction <- numeric(B_permutations)
for (b in seq_len(B_permutations)) {
  permuted_labels <- integer(nrow(doe_data))
  for (indices in clinic_indices) {
    permuted_labels[indices] <- sample(doe_data[["combination"]][indices])
  }
  permuted_interaction[b] <- interaction_statistic(
    permuted_labels, doe_data[["bp_reduction"]]
  )
}

omnibus_randomization_p <- (1 + sum(permuted_f >= observed_f)) /
  (B_permutations + 1)
interaction_randomization_p <-
  (1 + sum(abs(permuted_interaction) >= abs(observed_interaction))) /
  (B_permutations + 1)
randomization_result <- data.frame(
  Prespecified_statistic = c(
    "Interaction difference in differences (two-sided)",
    "Clinic-adjusted 3-df omnibus treatment F"
  ),
  Observed = c(observed_interaction, observed_f),
  Permutations = c(B_permutations, B_permutations),
  Fisher_sharp_null_p = c(interaction_randomization_p, omnibus_randomization_p)
)
knitr::kable(randomization_result, digits = 4,
  caption = "Two prespecified statistics under the same restricted randomization test of the global sharp null")

Two prespecified statistics under the same restricted randomization test of the global sharp null

Prespecified_statistic	Observed	Permutations	Fisher_sharp_null_p
Interaction difference in differences (two-sided)	2.25	999	0.158
Clinic-adjusted 3-df omnibus treatment F	25.73	999	0.001

The plus-one calculation avoids a zero Monte Carlo $p$-value. We permuted the four treatment labels jointly, not coaching and monitoring independently, and never moved a label across clinics. Both rows test the same sharp global null, but their sensitivity differs: the interaction statistic focuses on additive interaction, whereas the omnibus $F$ responds to any cell-mean difference and readily detects the large main effects here. Neither row is an exact test of an isolated weak or average-zero interaction null when main effects may exist.

old_par <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(
  permuted_interaction, breaks = 28, col = palette_doe["light"],
  border = "white", xlab = "Permuted interaction (mmHg)",
  main = "Interaction statistic"
)
abline(v = c(-abs(observed_interaction), abs(observed_interaction)),
       col = palette_doe["orange"], lwd = 2, lty = 2)
hist(
  permuted_f, breaks = 28, col = palette_doe["light"],
  border = "white", xlab = "Permuted omnibus F statistic",
  main = "Omnibus treatment F"
)
usr <- par("usr")
arrows(
  x0 = usr[2] - 0.04 * diff(usr[1:2]), y0 = usr[3] + 0.12 * diff(usr[3:4]),
  x1 = usr[2], y1 = usr[3] + 0.12 * diff(usr[3:4]),
  col = palette_doe["orange"], lwd = 2, length = 0.08
)
text(
  usr[2] - 0.04 * diff(usr[1:2]), usr[3] + 0.22 * diff(usr[3:4]),
  labels = sprintf("Observed F = %.2f (off scale)", observed_f),
  col = palette_doe["orange"], adj = 1, cex = 0.82
)

Restricted randomization distributions for the focused interaction contrast and the three-degree-of-freedom omnibus treatment F under the same global sharp null.

par(old_par)

8 Plan precision and power before collecting outcomes

8.1 What a useful calculation must specify

Power is a property of a proposed design under assumptions—not a quality score computed after seeing a $p$-value. A defensible plan names:

• one primary estimand and clinically meaningful target effect;
• anticipated outcome variance and its evidence source;
• factor allocation, blocks, correlation, attrition, and nonadherence;
• the intended model, contrast, significance level, multiplicity plan, and degrees of freedom;
• sensitivity to optimistic and pessimistic scenarios; and
• operational constraints such as clinic recruitment and minimum block size.

8.2 One-way ANOVA approximation

power.anova.test() treats the four combinations as independent groups with a common within-group variance. It can screen an omnibus four-cell question, but it ignores clinic blocking and does not target a particular factorial contrast. The rough calculation adds the prespecified clinic variance to the individual variance so it does not silently omit clinic-level outcome variation; it still misrepresents within-clinic correlation and the precision gain from blocking.

anova_power <- power.anova.test(
  groups = 4, n = NULL,
  between.var = var(unname(true_cell_means)),
  within.var = 4.5^2 + 2.2^2,
  sig.level = 0.05, power = 0.80
)
anova_power_summary <- data.frame(
  Target = "Four-group omnibus ANOVA approximation",
  Required_per_cell = ceiling(anova_power[["n"]]),
  Approximate_total = 4 * ceiling(anova_power[["n"]]),
  Assumed_within_SD = sqrt(4.5^2 + 2.2^2)
)
knitr::kable(anova_power_summary, digits = 2,
  caption = "Coarse four-group power approximation using prespecified simulation values")

Coarse four-group power approximation using prespecified simulation values

Target	Required_per_cell	Approximate_total	Assumed_within_SD
Four-group omnibus ANOVA approximation	10	40	5.01

Because the result does not enforce multiples of 12 clinics or participants per clinic-cell, it cannot be copied directly into this RCBD protocol. Design feasibility and the exact primary contrast must determine rounding.

8.3 Factorial simulation for the interaction

For a design-aligned calculation, simulate the proposed 12-clinic complete-block schedule, clinic and participant variation, the prespecified cell means, and the intended clinic + coaching * home_monitoring analysis. The target below is the 2.8-mmHg interaction. These assumptions are fixed in advance; this is not post hoc power based on the observed estimate.

simulate_interaction_power <- function(replicates_per_clinic_cell,
                                       simulations = 400) {
  simulated_schedule <- expand.grid(
    replicate = seq_len(replicates_per_clinic_cell),
    combination = seq_len(4),
    clinic_id = seq_len(n_clinics)
  )
  simulated_schedule[["coaching"]] <- factor(
    c(0, 1, 0, 1)[simulated_schedule[["combination"]]],
    levels = 0:1, labels = c("No", "Yes")
  )
  simulated_schedule[["monitoring"]] <- factor(
    c(0, 0, 1, 1)[simulated_schedule[["combination"]]],
    levels = 0:1, labels = c("No", "Yes")
  )
  simulated_schedule[["clinic"]] <- factor(simulated_schedule[["clinic_id"]])
  coach <- as.integer(simulated_schedule[["coaching"]] == "Yes")
  monitor <- as.integer(simulated_schedule[["monitoring"]] == "Yes")

  p_values <- numeric(simulations)
  for (s in seq_len(simulations)) {
    simulated_clinic_effect <- rnorm(n_clinics, 0, 2.2)
    simulated_outcome <- 4.0 + 3.0 * coach + 2.2 * monitor +
      2.8 * coach * monitor +
      simulated_clinic_effect[simulated_schedule[["clinic_id"]]] +
      rnorm(nrow(simulated_schedule), 0, 4.5)
    simulated_fit <- lm(
      simulated_outcome ~ clinic + coaching * monitoring,
      data = simulated_schedule
    )
    p_values[s] <- coef(summary(simulated_fit))[
      "coachingYes:monitoringYes", "Pr(>|t|)"
    ]
  }
  mean(p_values < 0.05)
}

set.seed(20260916)
replication_grid <- c(2, 4, 6, 8)
simulated_power <- vapply(
  replication_grid, simulate_interaction_power, numeric(1)
)
n_per_cell <- n_clinics * replication_grid
interaction_se_normal <- 2 * 4.5 / sqrt(n_per_cell)
noncentrality <- 2.8 / interaction_se_normal
z_critical <- qnorm(0.975)
normal_power <- pnorm(-z_critical - noncentrality) +
  1 - pnorm(z_critical - noncentrality)

factorial_power <- data.frame(
  Replicates_per_clinic_cell = replication_grid,
  Total_participants = 4 * n_per_cell,
  Simulated_power = simulated_power,
  Monte_Carlo_SE = sqrt(simulated_power * (1 - simulated_power) / 400),
  Normal_approximation = normal_power
)
knitr::kable(factorial_power, digits = 3,
  caption = "Interaction power under a 12-clinic balanced factorial design")

Interaction power under a 12-clinic balanced factorial design

Replicates_per_clinic_cell	Total_participants	Simulated_power	Monte_Carlo_SE	Normal_approximation
2	96	0.348	0.024	0.332
4	192	0.565	0.025	0.578
6	288	0.760	0.021	0.752
8	384	0.843	0.018	0.862

plot(
  factorial_power[["Total_participants"]], factorial_power[["Simulated_power"]],
  type = "b", pch = 19, lwd = 2, ylim = c(0, 1),
  col = palette_doe["teal"], xlab = "Total participants",
  ylab = "Power", main = "Power for the factorial interaction"
)
lines(
  factorial_power[["Total_participants"]], factorial_power[["Normal_approximation"]],
  type = "b", pch = 1, lwd = 2, col = palette_doe["orange"]
)
abline(h = 0.80, lty = 2, col = palette_doe["gray"])
legend(
  "bottomright", legend = c("Simulation", "Normal approximation", "80% target"),
  col = c(palette_doe["teal"], palette_doe["orange"], palette_doe["gray"]),
  pch = c(19, 1, NA), lty = c(1, 1, 2), bty = "n"
)

Simulated and normal-approximation power for the prespecified 2.8-mmHg interaction as total sample size increases in multiples allowed by 12 complete clinic blocks.

Simulation uncertainty remains: with 400 repetitions, a power estimate near 0.80 has Monte Carlo standard error about $\sqrt{0.8(0.2)/400}=0.02$. A protocol calculation should use more repetitions, vary all uncertain inputs, inflate recruitment for realistic outcome loss, and consider whether the interaction or a marginal main effect is truly primary.

9 Missingness, adherence, and the treatment-policy question

9.1 Preserve intention-to-treat

The randomized assignment is a baseline variable and must never be overwritten by treatment received. The primary intention-to-treat (ITT) comparison estimates the effect of assignment under real adherence in the trial. Per-protocol and as-treated analyses condition on postrandomization behavior and can be confounded.

set.seed(20260917)
adherence_example <- doe_data
adherence_example[["received_coaching"]] <- factor(
  ifelse(
    adherence_example[["coaching"]] == "Yes",
    rbinom(nrow(adherence_example), 1, 0.86),
    rbinom(nrow(adherence_example), 1, 0.08)
  ),
  0:1, c("No", "Yes")
)
adherence_table <- xtabs(~ coaching + received_coaching, data = adherence_example)
knitr::kable(adherence_table,
  caption = "Illustrative assigned versus received coaching; the ITT model retains assignment")

Illustrative assigned versus received coaching; the ITT model retains assignment

	No	Yes
No	89	7
Yes	6	90

Nonadherence is an outcome of the implementation process, not a reason to delete participants. If a per-protocol causal estimand is important, specify it in advance and use methods whose assumptions address postrandomization confounding, such as inverse-probability weighting or instrumental-variable approaches where justified.

Missing outcomes are different from nonadherence. Prevention begins with follow-up procedures and masked outcome collection. Report the number and reasons by randomized arm, define the missing-data estimand, and use sensitivity analyses for plausible departures from missing-at-random assumptions. Complete-case analysis can lose precision and be biased; simple last-observation-carried-forward or best/worst imputation does not generally solve the problem.

Ethics and design quality reinforce each other Clinical equipoise, a meaningful comparator, proportionate burden, informed consent, privacy protections, adverse-event monitoring, and a justified sample size are design requirements. Underpowered studies can expose participants without answering the question; excessive sample size can expose more participants than necessary. Subgroup and equity analyses should be planned with adequate representation and careful interpretation.

10 Avoid pseudoreplication and match analysis to assignment

10.1 The row count can be misleading

Suppose 12 clinics, rather than 192 participants, were randomized to a clinic-wide coaching policy. Outcomes from patients in the same clinic share one treatment assignment. Treating every patient as independently randomized would inflate the effective sample size: pseudoreplication. The treatment degrees of freedom and uncertainty must reflect 12 randomized clinics, along with intracluster correlation and small-cluster corrections where appropriate.

Technical replicates, repeated blood-pressure readings, bilateral organs, multiple time points, and multiple specimens can improve measurement but do not create new independent treatment assignments. Average or model subsamples at their proper level; never use them to manufacture experimental replication.

10.2 When another design is needed

Design	Assignment structure	Best use	Essential analysis boundary
Cluster-randomized	Groups receive treatments	Contamination, policy, or group delivery	Clusters are experimental units; plan ICC and enough clusters
Split-plot	Hard-to-change factor assigned to whole plots; another within plots	Multistage implementation or laboratory processing	Each effect uses its correct whole-plot or subplot error stratum
Crossover	Units receive sequences of treatments	Stable chronic conditions with reversible effects	Model period and within-person correlation; address carryover and washout
Latin square	Treatments balanced across two blocking dimensions	Two strong nuisance gradients	Basic square estimates treatment after rows/columns; interactions need replication
Repeated measures	Outcomes recur after one assignment	Trajectories and durability	Model within-unit correlation and time-by-treatment; visits are not new assignments
Incomplete block	Blocks cannot contain every treatment	Many treatments or capacity limits	Use the actual incidence structure and connected-design estimability

In crossover studies, dropout, period effects, and carryover can compromise a simple paired analysis. In split-plot studies, testing a whole-plot intervention against participant-level residual error is pseudoreplication. In Latin squares, a single unreplicated square cannot freely estimate all treatment-by-block interactions. Design diagrams and randomization records should accompany the analysis plan.

11 Outcomes beyond a continuous Gaussian response

The randomization principles remain, but the estimand and model scale change. A binary outcome may use a risk difference, risk ratio, or odds ratio; a count outcome may use a rate ratio with exposure time; repeated or clustered versions need appropriate correlation methods.

# Binary teaching outcome: at least an 8-mmHg reduction. Dichotomizing a
# continuous primary outcome loses information, so this would need a rationale.
outcome_example <- doe_data
outcome_example[["clinically_meaningful"]] <- outcome_example[["bp_reduction"]] >= 8
binary_fit <- glm(
  clinically_meaningful ~ clinic + coaching * home_monitoring,
  family = binomial, data = outcome_example
)

# Count teaching outcome with unequal observation time.
set.seed(20260918)
outcome_example[["followup_years"]] <- runif(nrow(outcome_example), 0.75, 1.00)
event_rate <- exp(-0.1 - 0.25 * (outcome_example[["coaching"]] == "Yes") -
                    0.15 * (outcome_example[["home_monitoring"]] == "Yes"))
outcome_example[["urgent_visits"]] <- rpois(
  nrow(outcome_example), outcome_example[["followup_years"]] * event_rate
)
count_fit <- glm(
  urgent_visits ~ clinic + coaching * home_monitoring +
    offset(log(followup_years)),
  family = poisson, data = outcome_example
)

outcome_model_examples <- data.frame(
  Outcome = c("Continuous BP reduction", "Binary meaningful reduction", "Urgent-visit count"),
  Model = c("Gaussian linear model", "Binomial logistic model", "Poisson log-rate model"),
  Natural_effect_scale = c("Mean difference", "Odds ratio in this model", "Rate ratio")
)
knitr::kable(outcome_model_examples,
  caption = "Outcome type changes the model and effect scale, not the randomization record")

Outcome type changes the model and effect scale, not the randomization record

Outcome	Model	Natural_effect_scale
Continuous BP reduction	Gaussian linear model	Mean difference
Binary meaningful reduction	Binomial logistic model	Odds ratio in this model
Urgent-visit count	Poisson log-rate model	Rate ratio

For common binary outcomes, odds ratios can be far from risk ratios; adjusted probabilities or risk differences may be easier to interpret. For counts, inspect overdispersion and use the log of person-time as an offset. For all generalized models, interaction is scale-specific, and model coefficients may not equal marginal causal contrasts. Standardize predicted outcomes across the planned target population when marginal effects are required.

12 Applications across public health and biomedicine

12.1 Translating the framework

Application	Factors and blocks	Primary design concern
Community hypertension program	Coaching × home monitoring; clinic blocks	Implementation fidelity and contamination
Vaccine formulation study	Antigen dose × adjuvant; site or production-lot blocks	Safety, dose-response, and multiplicity
Health-message experiment	Message frame × delivery channel; demographic strata	Interference and measurement of engagement
Laboratory assay optimization	Temperature × reagent concentration; plate/day blocks	Random run order and technical versus biological replication
Hospital implementation trial	Decision support × audit feedback; hospital clusters	Cluster-level assignment and few-cluster inference
Environmental intervention	Filtration × education; neighborhood blocks	Spillover, exposure measurement, and seasonal blocking

Factorial designs are especially attractive when interventions can be combined and resources must answer more than one question. They are inappropriate when combinations are unsafe, operationally impossible, or scientifically meaningless. Fractional factorial and response-surface designs can screen many factors or optimize doses, but aliasing, curvature, sequential experimentation, and confirmatory validation must be planned explicitly.

13 A reproducible analysis and reporting workflow

13.1 Selected results from the simulated trial

case_summary <- data.frame(
  Result = c(
    "Participants", "Clinics", "Participants per factorial cell",
    "Estimated marginal coaching effect (mmHg)",
    "Estimated marginal monitoring effect (mmHg)",
    "Estimated interaction (mmHg)",
    "Global sharp-null 3-df randomization p-value"
  ),
  Value = c(
    nrow(doe_data), nlevels(doe_data[["clinic"]]),
    min(as.vector(table(doe_data[["treatment_combination"]]))),
    contrast_results[["Estimate"]][1], contrast_results[["Estimate"]][2],
    contrast_results[["Estimate"]][3], omnibus_randomization_p
  )
)
knitr::kable(case_summary, digits = 3,
  caption = "Selected design and analysis results from the simulated experiment")

Selected design and analysis results from the simulated experiment

Result	Value
Participants	192.000
Clinics	12.000
Participants per factorial cell	48.000
Estimated marginal coaching effect (mmHg)	3.934
Estimated marginal monitoring effect (mmHg)	4.341
Estimated interaction (mmHg)	2.255
Global sharp-null 3-df randomization p-value	0.001

13.1.1 Auditable reporting template

We conducted a [CRD/RCBD/factorial/cluster/split-plot/crossover] experiment among [population, setting, and dates]. The experimental unit was [unit]; [factor levels] were assigned using [randomization procedure, restrictions, allocation ratio, and concealment], with [block/stratum] defined before assignment. The primary outcome was [definition, timing, direction], and the primary estimand was [explicit cell-mean contrast and population weighting]. We analyzed participants according to [assignment/other strategy] using [model], including [blocks, factors, interactions], and reported [effect scale] with 95% confidence intervals. Multiplicity was handled using [plan]. We assessed [missingness, adherence, residuals, dependence, influence, randomization-based sensitivity]. Limitations include [measurement, attrition, interference, generalizability, power, and design-specific constraints].

13.2 Analysis checklist

1. Lock the question, experimental unit, outcome time, estimand, and clinically meaningful effect.
2. Document eligibility, recruitment, intervention versions, control condition, and interference risks.
3. Choose factors, levels, allocation ratios, blocks, and any randomization restrictions.
4. Generate and secure a reproducible assignment schedule; preserve allocation concealment.
5. Plan blinding, measurement, adherence support, follow-up, harms monitoring, and data quality.
6. Calculate precision and power for the primary contrast under sensitivity scenarios.
7. Preregister primary and secondary outcomes, contrasts, exclusions, missing-data methods, and multiplicity.
8. Audit assignment and data integrity without testing baseline differences mechanically.
9. Analyze on the assignment scale using the design structure and correct experimental unit.
10. Report cell summaries, prespecified contrasts, uncertainty, diagnostics, attrition, deviations, and limits.

14 Common Errors at a Glance

Error	Why it is wrong	Better approach
Call every measured row an experimental unit	Measurements can share one assignment	Trace the level of independent randomization
Randomize globally but analyze as blocked	Analysis cannot invent restrictions not used	Record and reproduce the actual assignment mechanism
Shuffle coaching and monitoring separately in a randomization test	This may create assignments the design never allowed	Jointly permute four-arm labels within clinic
Interpret coachingYes as a marginal main effect	Under 0/1 coding it is simple at monitoring = No	Estimate an explicit four-cell contrast
Predict every cell at the reference clinic	The intercept encodes an arbitrary clinic	Average model-matrix rows across all target clinics
Drop the interaction after seeing its $p$-value	Data-dependent selection changes the estimand and uncertainty	Follow the prespecified hierarchy and report cell effects
Treat $p>0.05$ as proof of no effect	Imprecision is not equivalence	Report estimates, intervals, and meaningful bounds
Use repeated readings as independent replicates	They do not receive independent assignments	Model or summarize subsamples within experimental units
Exclude nonadherent participants from the primary analysis	Adherence is postrandomization	Preserve ITT and label secondary causal estimands
Calculate observed power	It adds no information beyond estimate and $p$-value	Plan prospective power across plausible scenarios
Ignore missing outcomes	Randomization does not eliminate attrition bias	Prevent loss, report flow, and perform sensitivity analyses
Generalize automatically from randomized participants	Randomization supports internal comparison, not sampling representativeness	Define the target population and assess transportability

15 Exercises and Answers

1. In this trial, are clinics, participants, or BP readings the experimental units? Why?
2. Under 0/1 coding, what does coachingYes estimate when an interaction is in the model?
3. Using cell means $(4.0,7.0,6.2,12.0)$, calculate the marginal coaching and interaction effects.
4. Why should cell means be averaged over all 12 clinic design rows instead of evaluated at Clinic 01?
5. How must treatment labels be permuted for a valid randomization test of the global sharp null?
6. A 95% CI for an effect is $-0.4$ to $3.2$ mmHg. What can and cannot be concluded?
7. If coaching is assigned once per clinic, can 100 patients per clinic be treated as 100 independent coaching assignments?
8. Why is power.anova.test() only a rough approximation for the primary interaction here?

Show exercise answers

1. Participants are the experimental units because treatment combinations are independently assigned to participants within clinics. Clinics are blocks, and repeat BP readings would be observational subsamples unless separately randomized.
2. It estimates coaching’s simple effect when home monitoring is at its reference level, “No.” The marginal coaching effect is $\beta_A+\beta_{AB}/2$ under equal 50/50 averaging over monitoring.
3. Marginal coaching is $[(-4+7-6.2+12)/2]=4.4$ mmHg. Interaction is $4-7-6.2+12=2.8$ mmHg.
4. Clinic indicators shift expected outcomes. Evaluating at Clinic 01 gives a reference-clinic mean; averaging identical treatment profiles over all clinics gives the same clinic distribution to every cell and targets the study’s clinic mix.
5. Treat each four-level combination label as one object, shuffle labels only among participants in the same clinic, and retain four assignments per combination per clinic. Do not shuffle factor columns independently or across clinics.
6. Values from a small adverse effect through a potentially important benefit are compatible with the model and data. The interval does not prove no effect, equivalence, or a clinically important benefit; those require prespecified margins and adequate precision.
7. No. The clinic is the experimental unit for coaching; patients share its assignment. The independent treatment information comes from randomized clinics, with within-clinic outcomes correlated.
8. It treats four independent groups, ignores clinic blocking, and tests omnibus between-cell variation rather than the named interaction contrast. A design-aligned simulation better represents replication, blocks, analysis, and attrition.

16 Quick Reference

16.1 Core formulas and R entry points

Target	Formula or code	Interpretation
RCBD factorial model	lm(y ~ block + A * B)	Block-adjusted cells, main terms, and interaction
Coaching marginal effect	$(-\mu_{00}+\mu_{10}-\mu_{01}+\mu_{11})/2$	Coaching averaged equally over monitoring
Monitoring marginal effect	$(-\mu_{00}-\mu_{10}+\mu_{01}+\mu_{11})/2$	Monitoring averaged equally over coaching
Interaction	$\mu_{00}-\mu_{10}-\mu_{01}+\mu_{11}$	Difference between two simple effects
Linear contrast variance	$\operatorname{Var}(L\hat\beta)=L\widehat V L^T$	Uses full coefficient covariance
Pairwise cell comparisons	TukeyHSD(aov_fit)	Simultaneous pairwise intervals in a suitable ANOVA
Global randomization test	Permute joint arm labels within block	Tests the sharp null under the actual restriction
Four-group screening power	power.anova.test(...)	Omnibus approximation, not a factorial contrast
Binary outcome	glm(..., family = binomial)	Conditional odds on the logit scale
Count rate	glm(... + offset(log(time)), family = poisson)	Conditional event-rate comparison

16.2 Final design checklist

• Is the experimental unit defined from the randomization operation?
• Are factor levels, intervention versions, comparator, outcome timing, and estimand explicit?
• Does allocation concealment prevent foreknowledge, and is blinding used where feasible?
• Do blocks use preassignment information and contain the comparisons they are meant to improve?
• Is replication independent rather than technical or longitudinal pseudoreplication?
• Are the primary cell-mean contrast, effect scale, multiplicity, and missing-data plan prespecified?
• Does sample-size planning match the actual factorial, blocking, clustering, and attrition structure?
• Does the analysis retain assignment, blocks, hierarchy, and the correct experimental-unit error?
• Are marginal effects distinguished from reference-level simple effects?
• Are effect estimates, confidence intervals, absolute cell means, diagnostics, adherence, and attrition reported?
• Are interference, measurement, ethics, generalizability, and design deviations discussed?

Next directions for study

Next topics include covariate-adjusted randomized-trial analysis, ANCOVA with baseline outcomes, response-surface methodology, fractional factorial screening, optimal designs, sequential and adaptive designs, noninferiority and equivalence, group-sequential monitoring, stepped-wedge cluster trials, randomization-based confidence intervals, treatment-effect heterogeneity, mediation, principal stratification, and transportability.

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.60       cachem_1.1.0   
##  [6] knitr_1.51      htmltools_0.5.9 rmarkdown_2.31  lifecycle_1.0.5 cli_3.6.6      
## [11] sass_0.4.10     jquerylib_0.1.4 compiler_4.6.1  tools_4.6.1     evaluate_1.0.5 
## [16] bslib_0.12.0    yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

For teaching only: This module introduces experimental-design principles, factorial estimation, model and randomization-based inference, power, diagnostics, and reporting. It does not replace a protocol-specific statistical analysis plan, ethics review, trial registration, data-monitoring plan, or consultation with a biostatistician. All data and numerical results are simulated.