Randomized Complete Block Design (RCBD)

Controlling Heterogeneity to Increase Precision

Dr. Samuel B Fernandes

2026-01-01

Learning Objectives

By the end of this lecture, you will be able to:

Explain why and when to use blocking to control heterogeneity
Design and randomize a Randomized Complete Block Design (RCBD)
Analyze RCBD models
Interpret block and treatment effects from ANOVA tables

R packages used in this lecture

library(tidyverse) # Data manipulation and visualization
library(car) # ANOVA Type III sums of squares & regression diagnostics
library(emmeans) # Estimated marginal means & post-hoc comparisons (Tukey)
library(agridat) # Agricultural datasets (yates.oats, etc.)
library(easystats) # Comprehensive model diagnostics, normality, homogeneity
library(asbio) # Tukey's additivity test for checking RCBD assumptions

Part 1: Why Block?

The Challenge: Hidden Variation

The Problem: Even within a single experiment, experimental units are not identical.

Common sources of heterogeneity in agricultural research:

Field trials: Soil texture gradients, drainage patterns, fertility zones, edge effects
Greenhouse: Bench position (light, temperature), proximity to doors/vents
Orchards: Tree age, row position, rootstock variation
Animal studies: Litter, age, initial weight, genetic background
Time effects: Planting date, harvest day, processing batch, weather conditions

Visual: The Problem with CRD

Code

set.seed(2026)
field_gradient <- data.frame(
    plot = 1:20,
    position = 1:20,
    treatment = sample(rep(c("A", "B", "C", "D"), each = 5)),
    soil_quality = seq(10, 5, length.out = 20) # Gradient from good to poor
) |>
    mutate(
        yield_base = 100 + soil_quality * 2,
        treatment_effect = case_when(
            treatment == "A" ~ 0,
            treatment == "B" ~ 3,
            treatment == "C" ~ 6,
            TRUE ~ 10
        ),
        yield = yield_base + treatment_effect + rnorm(20, 0, 3)
    )

ggplot(field_gradient, aes(x = position, y = yield, color = treatment, shape = treatment)) +
    geom_point(size = 4) +
    geom_smooth(aes(group = 1), method = "lm", se = FALSE, color = "gray50", linetype = "dashed") +
    labs(
        title = "CRD: Treatments randomized across field with soil gradient",
        subtitle = "Gray line shows underlying soil quality gradient (confounded with treatments)",
        x = "Field Position", y = "Yield (kg/ha)", color = "Treatment", shape = "Treatment"
    ) +
    theme_minimal(base_size = 14) +
    theme(legend.position = "bottom")

Figure 1: CRD with hidden gradient: treatments randomly assigned across heterogeneous field

Issue: Treatments at different positions experience different soil conditions → inflated error variance

The Solution: Blocking

Key Idea: Group similar experimental units into blocks; randomize treatments within each block.

Benefits of Blocking:

Removes known sources of variation from error term
Smaller MS_Error → more powerful F-tests
Fair comparison of treatments (each block gets all treatments)
Increased precision without increasing sample size

When to Block:

✅ Block when you can identify and group similar experimental units before the experiment

✅ Block on factors you’re not interested in testing (nuisance factors)

❌ Don’t block on factors you want to compare (those are treatments!)

Visual: RCBD Solution

Code

set.seed(2026)
field_rcbd <- data.frame(
    block = rep(1:5, each = 4),
    position = rep(c(1, 2, 3, 4), times = 5)
) |>
    group_by(block) |>
    mutate(treatment = sample(c("A", "B", "C", "D"))) |>
    ungroup() |>
    mutate(
        field_position = (block - 1) * 4 + position,
        soil_quality = seq(10, 5, length.out = 20)[field_position],
        yield_base = 100 + soil_quality * 2,
        treatment_effect = case_when(
            treatment == "A" ~ 0,
            treatment == "B" ~ 3,
            treatment == "C" ~ 6,
            TRUE ~ 10
        ),
        yield = yield_base + treatment_effect + rnorm(20, 0, 3),
        block = factor(block)
    )

ggplot(field_rcbd, aes(x = field_position, y = yield, color = treatment, shape = treatment)) +
    geom_point(size = 4) +
    geom_vline(xintercept = c(4.5, 8.5, 12.5, 16.5), linetype = "dotted", color = "black", linewidth = 0.8) +
    annotate("text",
        x = c(2.5, 6.5, 10.5, 14.5, 18.5), y = 129,
        label = paste("Block", 1:5), color = "black", size = 4
    ) +
    labs(
        title = "RCBD: Treatments randomized within blocks",
        subtitle = "Each block contains all 4 treatments; blocks account for soil gradient",
        x = "Field Position", y = "Yield (kg/ha)", color = "Treatment", shape = "Treatment"
    ) +
    theme_minimal(base_size = 14) +
    theme(legend.position = "bottom")

Figure 2: RCBD: Field divided into blocks; treatments randomized within each block

Result: Block-to-block variation removed from error term; treatments compared fairly within blocks

RCBD Model & Structure

The RCBD Statistical Model

\[ y_{ij} = \mu + \beta_j + \tau_i + \epsilon_{ij} \]

Where:

\(y_{ij}\) = response for treatment \(i\) in block \(j\)
\(\mu\) = overall mean
\(\beta_j\) = effect of block \(j\) (removes known heterogeneity)
\(\tau_i\) = effect of treatment \(i\) (what we want to test!)
\(\epsilon_{ij} \sim N(0, \sigma^2)\) = random error (within-block variability)

Key assumption: Additivity — block and treatment effects are additive (no interaction)

Additivity means: The effect of a treatment is the same across all blocks. Block effects shift all treatments up or down equally.

ANOVA Table: CRD vs RCBD

CRD (One-Way ANOVA)

Source	df	SS	MS	F
Treatment	\(t-1\)	SS_Trt	MS_Trt	F_Trt
Error	\(t(r-1)\)	SS_E	MS_E
Total	\(tr-1\)	SS_Total

Error includes:

Within-treatment variation
+ Block-to-block variation ❌

RCBD

Source	df	SS	MS	F
Block	\(b-1\)	SS_Block	MS_Block	F_Block
Treatment	\(t-1\)	SS_Trt	MS_Trt	F_Trt
Error	\((b-1)(t-1)\)	SS_E	MS_E
Total	\(bt-1\)	SS_Total

Error includes:

Within-block variation only ✅
Block variation removed!

Result: \(MS_{Error(RCBD)} < MS_{Error(CRD)}\) → more powerful tests

Assumptions for RCBD

Same as CRD, plus one more:

Normality: Errors \(\epsilon_{ij}\) are normally distributed
Homogeneity of variance: Equal variance across treatment-block combinations
Independence: Observations are independent (proper randomization)
Additivity (NEW!): No block × treatment interaction

Caution

Important: If interaction exists (non-additivity), we need a mixed model or analyze blocks separately! Or multiple reps/block to test interaction.

we can check it with asbio::tukey.add.test()!

Part 2: Examples & Applications

Example 1: Flower Vase Experiment

Research Question: Does the liquid in a vase affect how long cut flowers last?

Experimental Design:

Treatments (4): Tap water, Sugar water, Carbonated water, 7-Up
Blocks (3): Flower species (Rose, Daisy, Tulip)
Response: Days until flower wilts
Why block? Different species have different baseline vase life

Each species (block) gets all 4 liquid treatments → complete block design

Flower Vase: Design & Randomization

By hand or using the FielDHub Shiny app.

Code

set.seed(2026)
trt_levels <- c("Tap", "Sugar", "Carbonated", "7-Up")
b1 <- sample(trt_levels, 4)
b2 <- sample(trt_levels, 4)
b3 <- sample(trt_levels, 4)
block <- factor(rep(c("daisy", "rose", "tulip"), each = 4))
flower_plan <- data.frame(TypeFlower = block, FlowerNumber = 1:12, treatment = c(b1, b2, b3))
# write_csv(flower_plan, file = "RCBPlan.csv")

Randomization table shows: Each block (species) has all 4 treatments in random order

Flower Vase: Simulated Data & Analysis

Code

#|
# Simulate realistic data
flower_data <- flower_plan |>
    mutate(
        # Add species labels based on block
        species = c("Rose", "Daisy", "Tulip")[block],
        # Species baseline effect
        species_effect = case_when(
            species == "Rose" ~ 3,
            species == "Daisy" ~ -1,
            species == "Tulip" ~ 0,
            TRUE ~ 0
        ),
        # Treatment effect
        trt_effect = case_when(
            treatment == "Tap" ~ 0,
            treatment == "Sugar" ~ 3,
            treatment == "Carbonated" ~ 1,
            treatment == "7-Up" ~ 4,
            TRUE ~ 0
        ),
        days_wilted = 7 + species_effect + trt_effect + rnorm(n(), 0, 1.5)
    )

# Visualize: Boxplot by treatment, colored by block
ggplot(flower_data, aes(x = treatment, y = days_wilted, fill = species)) +
    geom_boxplot(alpha = 0.7) +
    geom_point(position = position_jitterdodge(jitter.width = 0.1), size = 2) +
    labs(
        title = "Flower Vase Life by Treatment and Species",
        x = "Vase Liquid", y = "Days Until Wilted", fill = "Flower Species"
    ) +
    theme_minimal(base_size = 14) +
    theme(legend.position = "bottom")

Flower Vase: ANOVA with Type III SS

Critical: Set contrasts before fitting the model for Type III SS!

Code

# Set sum-to-zero contrasts for Type III SS
options(contrasts = c("contr.sum", "contr.poly"))

# Fit RCBD model: y ~ block + treatment
model_flower <- lm(days_wilted ~ block + treatment, data = flower_data)

# Type III ANOVA (order-independent)
Anova(model_flower, type = "III")

#> Anova Table (Type III tests)
#> 
#> Response: days_wilted
#>             Sum Sq Df  F value    Pr(>F)    
#> (Intercept) 814.85  1 399.7067 1.016e-06 ***
#> block        29.50  2   7.2361   0.02517 *  
#> treatment    54.69  3   8.9421   0.01241 *  
#> Residuals    12.23  6                       
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretation:

Block (species) effect: p < 0.05 → species differ in baseline vase life
Treatment effect: p < 0.05 → liquid type affects vase life
Use smaller MS_Error from RCBD (not from ignoring blocks!)

Flower Vase: Model Diagnostics with easystats

Code

# Check all assumptions at once
check_model(model_flower)

easystats interpretation:

check_model() produces 6 diagnostic plots automatically
Look for: normality (Q-Q plot), homogeneity (scale-location), outliers, influential points

Flower Vase: Post-hoc Comparisons

Code

# Estimated marginal means by treatment
emm_flower <- emmeans(model_flower, ~treatment)

# Pairwise comparisons with Tukey adjustment
pairs(emm_flower, adjust = "tukey")

#>  contrast            estimate   SE df t.ratio p.value
#>  (7-Up) - Carbonated     5.86 1.17  6   5.025  0.0095
#>  (7-Up) - Sugar          2.19 1.17  6   1.877  0.3281
#>  (7-Up) - Tap            3.65 1.17  6   3.132  0.0735
#>  Carbonated - Sugar     -3.67 1.17  6  -3.148  0.0722
#>  Carbonated - Tap       -2.21 1.17  6  -1.893  0.3222
#>  Sugar - Tap             1.46 1.17  6   1.255  0.6191
#> 
#> Results are averaged over the levels of: block 
#> P value adjustment: tukey method for comparing a family of 4 estimates

Code

# Visualize with confidence intervals
plot(emm_flower, comparisons = TRUE) +
    labs(title = "Treatment Means with 95% CIs") +
    theme_minimal(base_size = 14)

Results: 7-Up and Sugar water significantly extend vase life vs. Tap water

Check Additivity with Interaction Plot

Code

# Interaction plot to check additivity assumption
ggplot(flower_data, aes(x = treatment, y = days_wilted, color = species, group = species)) +
    stat_summary(fun = mean, geom = "point", size = 3) +
    stat_summary(fun = mean, geom = "line", linewidth = 1) +
    labs(
        title = "Interaction Plot: Check for Additivity",
        subtitle = "Parallel lines indicate no block × treatment interaction (good!)",
        x = "Treatment", y = "Mean Days Until Wilted", color = "Species"
    ) +
    theme_minimal(base_size = 14) +
    theme(legend.position = "bottom")

Interpretation: Lines are roughly parallel → additivity assumption satisfied ?

Check Additivity with Tukey’s Additivity Test

tukey.add.test(flower_data$days_wilted, flower_data$treatment, flower_data$species)

#> 
#> Tukey's one df test for additivity 
#> F = 0.2870526   Denom df = 5    p-value = 0.6150809

Example 3: Real Agricultural Data

Example 3: Oats Yield Trial (agridat)

Using the classic Yates oats dataset from agridat package:

Code

# agridat package already loaded in setup
data(yates.oats)

# Examine structure
str(yates.oats)

#> 'data.frame':    72 obs. of  8 variables:
#>  $ row  : int  16 12 3 14 8 5 15 11 3 14 ...
#>  $ col  : int  3 4 3 1 2 2 4 4 4 2 ...
#>  $ yield: int  80 60 89 117 64 70 82 102 82 114 ...
#>  $ nitro: num  0 0 0 0 0 0 0.2 0.2 0.2 0.2 ...
#>  $ gen  : Factor w/ 3 levels "GoldenRain","Marvellous",..: 1 1 1 1 1 1 1 1 1 1 ...
#>  $ block: Factor w/ 6 levels "B1","B2","B3",..: 1 2 3 4 5 6 1 2 3 4 ...
#>  $ grain: num  20 15 22.2 29.2 16 ...
#>  $ straw: num  28 25 40.5 28.8 32 ...

Code

# Summary
head(yates.oats)

#>   row col yield nitro        gen block grain straw
#> 1  16   3    80     0 GoldenRain    B1 20.00 28.00
#> 2  12   4    60     0 GoldenRain    B2 15.00 25.00
#> 3   3   3    89     0 GoldenRain    B3 22.25 40.50
#> 4  14   1   117     0 GoldenRain    B4 29.25 28.75
#> 5   8   2    64     0 GoldenRain    B5 16.00 32.00
#> 6   5   2    70     0 GoldenRain    B6 17.50 27.25

Study details:

Treatments: Nitrogen levels (0, 0.2, 0.4, 0.6 cwt/acre)
Blocks: 6 blocks
Varieties: 3 oat varieties (Golden Rain, Marvellous, Victory)
Response: Grain yield

Oats: Visualize Block & Treatment Structure

Code

# Focus on nitrogen effect within one variety
oats_subset <- yates.oats

ggplot(oats_subset, aes(x = factor(nitro), y = yield, fill = block)) +
    geom_boxplot(alpha = 0.7) +
    geom_point(position = position_jitterdodge(jitter.width = 0.1), size = 2) +
    labs(
        title = "Oats Yield by Nitrogen Level",
        x = "Nitrogen (cwt/acre)", y = "Yield (bushels/acre)", fill = "Block"
    ) +
    theme_minimal(base_size = 14) +
    theme(legend.position = "bottom")

Oats: RCBD Analysis

Code

# Set contrasts
options(contrasts = c("contr.sum", "contr.poly"))

# Fit RCBD model (include variety as additional factor)
model_oats <- lm(yield ~ block + gen + factor(nitro), data = oats_subset)

# Type III ANOVA
Anova(model_oats, type = "III")

#> Anova Table (Type III tests)
#> 
#> Response: yield
#>               Sum Sq Df   F value    Pr(>F)    
#> (Intercept)   778336  1 3319.2914 < 2.2e-16 ***
#> block          15875  5   13.5403 6.906e-09 ***
#> gen             1786  2    3.8091   0.02762 *  
#> factor(nitro)  20020  3   28.4598 1.239e-11 ***
#> Residuals      14304 61                        
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# Post-hoc for nitrogen effect
emm_oats <- emmeans(model_oats, ~nitro)
pairs(emm_oats, adjust = "tukey")

#>  contrast            estimate  SE df t.ratio p.value
#>  nitro0 - nitro0.2     -19.50 5.1 61  -3.820  0.0017
#>  nitro0 - nitro0.4     -34.83 5.1 61  -6.824  <.0001
#>  nitro0 - nitro0.6     -44.00 5.1 61  -8.620  <.0001
#>  nitro0.2 - nitro0.4   -15.33 5.1 61  -3.004  0.0197
#>  nitro0.2 - nitro0.6   -24.50 5.1 61  -4.800  0.0001
#>  nitro0.4 - nitro0.6    -9.17 5.1 61  -1.796  0.2852
#> 
#> Results are averaged over the levels of: block, gen 
#> P value adjustment: tukey method for comparing a family of 4 estimates

Result: Nitrogen significantly increases oats yield; blocking on field position reduced error

Oats: Model Performance Report

Code

# easystats package already loaded in setup; use report() function
report(model_oats)

#> We fitted a linear model (estimated using OLS) to predict yield with block, gen
#> and nitro (formula: yield ~ block + gen + factor(nitro)). The model explains a
#> statistically significant and substantial proportion of variance (R2 = 0.72,
#> F(10, 61) = 16.07, p < .001, adj. R2 = 0.68). The model's intercept,
#> corresponding to block = , gen = and nitro = 0, is at 103.97 (95% CI [100.36,
#> 107.58], t(61) = 57.61, p < .001). Within this model:
#> 
#>   - The effect of block [1] is statistically significant and negative (beta =
#> -13.06, 95% CI [-21.12, -4.99], t(61) = -3.24, p = 0.002; Std. beta = -0.48,
#> 95% CI [-0.78, -0.18])
#>   - The effect of block [2] is statistically non-significant and negative (beta =
#> -8.06, 95% CI [-16.12, 0.01], t(61) = -2.00, p = 0.050; Std. beta = -0.30, 95%
#> CI [-0.60, 5.02e-04])
#>   - The effect of block [3] is statistically non-significant and negative (beta =
#> -7.72, 95% CI [-15.79, 0.35], t(61) = -1.91, p = 0.060; Std. beta = -0.29, 95%
#> CI [-0.58, 0.01])
#>   - The effect of block [4] is statistically significant and positive (beta =
#> 31.36, 95% CI [23.29, 39.43], t(61) = 7.77, p < .001; Std. beta = 1.16, 95% CI
#> [0.86, 1.46])
#>   - The effect of block [5] is statistically non-significant and negative (beta =
#> -5.81, 95% CI [-13.87, 2.26], t(61) = -1.44, p = 0.155; Std. beta = -0.21, 95%
#> CI [-0.51, 0.08])
#>   - The effect of gen [1] is statistically non-significant and positive (beta =
#> 0.53, 95% CI [-4.58, 5.63], t(61) = 0.21, p = 0.837; Std. beta = 0.02, 95% CI
#> [-0.17, 0.21])
#>   - The effect of gen [2] is statistically significant and positive (beta = 5.82,
#> 95% CI [0.72, 10.92], t(61) = 2.28, p = 0.026; Std. beta = 0.22, 95% CI [0.03,
#> 0.40])
#>   - The effect of nitro [1] is statistically significant and negative (beta =
#> -24.58, 95% CI [-30.83, -18.33], t(61) = -7.86, p < .001; Std. beta = -0.91,
#> 95% CI [-1.14, -0.68])
#>   - The effect of nitro [2] is statistically non-significant and negative (beta =
#> -5.08, 95% CI [-11.33, 1.17], t(61) = -1.63, p = 0.109; Std. beta = -0.19, 95%
#> CI [-0.42, 0.04])
#>   - The effect of nitro [3] is statistically significant and positive (beta =
#> 10.25, 95% CI [4.00, 16.50], t(61) = 3.28, p = 0.002; Std. beta = 0.38, 95% CI
#> [0.15, 0.61])
#> 
#> Standardized parameters were obtained by fitting the model on a standardized
#> version of the dataset. 95% Confidence Intervals (CIs) and p-values were
#> computed using a Wald t-distribution approximation.

The report() function generates a publication-ready summary in plain English!

Example 4: Food Science—Noodle Quality Evaluation

Example 4: Sensory Panel Study of Noodle Varieties

Research question: Do cooking time (1, 2, 3, 4 minutes) affect sensory quality (texture, firmness) across different noodle varieties?

Why RCBD?

Panelists are “blocks” — each panelist has inherent sensory sensitivity differences that shouldn’t confound treatment effects
Treatments: Cooking times (fixed, of interest)
Blocks: Panelist judges (random, control for individual variation)

Noodle Study: Design & Data

Code

# Simulate a sensory panel study: 12 judges rating noodles at 4 cooking times
set.seed(2025)

noodle_data <- expand_grid(
    judge = 1:12,
    cooking_time = c(1, 2, 3, 4) # minutes
) |>
    mutate(
        # Judge baseline effect (consistent across treatments)
        judge_bias = rep(c(-2, -1, 0, 0.5, 1, 1.5, -1.5, 0, 0.5, 1, -0.5, 0.5), each = 4),

        # Cooking time effect on texture/firmness score (1-9 scale)
        time_effect = case_when(
            cooking_time == 1 ~ -1.5, # Too firm
            cooking_time == 2 ~ 2.5, # Optimal
            cooking_time == 3 ~ 1.0, # Slightly soft
            cooking_time == 4 ~ -2.0, # Too soft
            TRUE ~ 0
        ),

        # Texture firmness score (1-9 scale)
        firmness = 5 + judge_bias + time_effect + rnorm(n(), 0, 0.8)
    ) |>
    select(-judge_bias, -time_effect)

head(noodle_data, 12)

#> # A tibble: 12 × 3
#>    judge cooking_time firmness
#>    <int>        <dbl>    <dbl>
#>  1     1            1     2.00
#>  2     1            2     5.53
#>  3     1            3     4.62
#>  4     1            4     2.02
#>  5     2            1     2.80
#>  6     2            2     6.37
#>  7     2            3     5.32
#>  8     2            4     1.94
#>  9     3            1     3.22
#> 10     3            2     8.06
#> 11     3            3     5.68
#> 12     3            4     1.60

Noodle Study: Visualization

Code

ggplot(noodle_data, aes(x = factor(cooking_time), y = firmness)) +
    geom_boxplot(alpha = 0.6, fill = "steelblue") +
    geom_point(aes(color = factor(judge)), alpha = 0.5, size = 2) +
    geom_line(aes(group = judge, color = factor(judge)), alpha = 0.3) +
    labs(
        title = "Noodle Firmness Scores by Cooking Time (Panelists as Blocks)",
        x = "Cooking Time (minutes)",
        y = "Firmness Score (1-9)",
        color = "Judge ID"
    ) +
    theme_minimal(base_size = 13) +
    theme(legend.position = "right", panel.grid.major.x = element_blank())

Observation: Judges differ in baseline firmness ratings (blocking effect) but similar pattern across cooking times.

Noodle Study: RCBD Analysis

Code

# Set sum-to-zero contrasts BEFORE fitting model
options(contrasts = c("contr.sum", "contr.poly"))

# RCBD model: firmness ~ judge (block) + cooking_time (treatment)
model_noodle <- lm(firmness ~ factor(judge) + factor(cooking_time), data = noodle_data)

# Type III ANOVA
Anova(model_noodle, type = "III")

#> Anova Table (Type III tests)
#> 
#> Response: firmness
#>                       Sum Sq Df   F value    Pr(>F)    
#> (Intercept)          1249.77  1 1805.1085 < 2.2e-16 ***
#> factor(judge)          63.17 11    8.2951 9.877e-07 ***
#> factor(cooking_time)  171.17  3   82.4115 2.073e-15 ***
#> Residuals              22.85 33                        
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretation:

Judge effect (p < 0.001): Panelists have significantly different baseline firmness expectations → blocking was necessary
Cooking time effect (p < 0.001): Cooking duration significantly affects noodle firmness → 2 minutes is optimal
Compare: If we’d ignored judges (CRD), judge variation would inflate error and reduce power

Noodle Study: Post-hoc Comparisons

Code

# Estimated marginal means for cooking times
emm_noodle <- emmeans(model_noodle, ~ factor(cooking_time))

# Tukey-adjusted pairwise differences
pairs(emm_noodle, adjust = "tukey")

#>  contrast                      estimate   SE df t.ratio p.value
#>  cooking_time1 - cooking_time2   -4.362 0.34 33 -12.840  <.0001
#>  cooking_time1 - cooking_time3   -2.665 0.34 33  -7.845  <.0001
#>  cooking_time1 - cooking_time4    0.133 0.34 33   0.393  0.9791
#>  cooking_time2 - cooking_time3    1.697 0.34 33   4.995  0.0001
#>  cooking_time2 - cooking_time4    4.495 0.34 33  13.233  <.0001
#>  cooking_time3 - cooking_time4    2.798 0.34 33   8.238  <.0001
#> 
#> Results are averaged over the levels of: judge 
#> P value adjustment: tukey method for comparing a family of 4 estimates

Code

# Plot with confidence intervals
plot(emm_noodle, comparisons = TRUE) +
    labs(
        title = "Noodle Firmness by Cooking Time (95% CI)",
        x = "Estimated Firmness Score"
    ) +
    theme_minimal(base_size = 13)

Results:

2 min vs. 1 min: Significantly higher firmness (p < 0.01) — cooking time improves texture
2 min vs. 3 min: Slightly higher firmness (p = 0.08) — approaching over-softness
2 min vs. 4 min: Much higher firmness (p < 0.001) — clear over-cooking damage

Practical Decision: Recommend 2-minute cooking time for optimal sensory response across all panelist groups.

Connection to Factorial Designs

RCBD as Foundation for Factorial Designs

Key Connection: RCBD with one blocking factor and one treatment factor is a two-way design.

In Week 5, we’ll extend this to:

Factorial CRD: Two or more treatment factors, no blocking
Factorial RCBD: Two or more treatment factors + blocking
Interactions: When treatment effects depend on levels of another factor

Example: Herbicide × Application Timing factorial in RCBD blocks

Preview: Two-Way ANOVA Structure

Model for Factorial RCBD:

\[ y_{ijk} = \mu + \beta_k + \alpha_i + \gamma_j + (\alpha\gamma)_{ij} + \epsilon_{ijk} \]

Where:

\(\beta_k\) = block effect
\(\alpha_i\) = effect of factor A (e.g., herbicide)
\(\gamma_j\) = effect of factor B (e.g., application timing)
\((\alpha\gamma)_{ij}\) = interaction between A and B

Today’s RCBD: Special case where \(\alpha\) = treatment, \(\beta\) = block, no interaction term

Summary & Key Takeaways

Key Concepts to Remember

Blocking removes known heterogeneity from the error term → smaller MS_Error → more powerful tests
RCBD requires complete blocks: Each block gets all treatments
Randomize treatments within blocks, not across the entire experiment
Set contrasts before fitting: options(contrasts = c("contr.sum", "contr.poly")) for Type III SS
Check additivity: Use interaction plots to verify no block × treatment interaction
Use easystats for diagnostics: performance::check_model() provides comprehensive assumption checks
RCBD extends naturally to factorial designs with multiple treatment factors

Common Mistakes to Avoid

❌ Don’t do this:

Forget to include block in the model (analyzing as CRD)
Block on factors you want to test (treatments should not be blocks!)
Use Type I SS without setting contrasts (order-dependent results)
Ignore additivity assumption (check interaction plots!)
Analyze with unbalanced blocks without checking missing data patterns

✅ Best practices:

Always visualize data by block and treatment before analysis
Use Anova(..., type = "III") from car package
Check assumptions with check_model() from performance
Report both F-tests and post-hoc comparisons with adjusted p-values

References & Resources

Textbook:

Oehlert (2010), A First Course in Design and Analysis of Experiments

R Packages:

agridat: 100+ agricultural datasets with RCBD examples
car: Companion to Applied Regression (Type III ANOVA)
emmeans: Estimated marginal means (Lenth, 2022)
easystats: Collection of packages for intuitive model analysis