Latin squares and unequal group sizes

After the break

Load the packages for this workshop

Use the same modelling tools as before the break. readr reads the CSV files, dplyr prepares and counts records, and tidyr changes layouts and creates combinations. tibble creates the small simulated dataset used later.

library("lme4")
library("performance")
library("marginaleffects", warn.conflicts = FALSE)
library("ggplot2")
library("readr")
library("dplyr")
library("tidyr")
library("tibble")

Compare four feeds for broilers

Sixteen pens occupy four rows and four columns in a broiler house. Each pen holds the same number of birds of the same age and receives one of four feeds, A–D, throughout the trial.

Light levels differ between rows; airflow differs between columns. We want to compare feeds while allowing for both sources of variation.

A Latin square places every feed once in each row and once in each column.

Read the pen layout

Each record gives one pen’s position, feed and mean daily gain. We record one mean daily gain per pen, in grams per bird per day, over the same trial period.

l <- read_csv("broiler-latin-square.csv", show_col_types = FALSE) |>
  mutate(
    row = factor(row, levels = paste0("Row", 1:4)),
    column = factor(column, levels = paste0("Col", 1:4)),
    feed = factor(feed, levels = c("A", "B", "C", "D"))
  )
l |>
  select(row, column, feed) |>
  pivot_wider(names_from = column, values_from = feed)

# A tibble: 4 × 5
  row   Col1  Col2  Col3  Col4 
  <fct> <fct> <fct> <fct> <fct>
1 Row1  B     D     A     C    
2 Row2  C     A     B     D    
3 Row3  A     C     D     B    
4 Row4  D     B     C     A    

Find all four pens receiving feed A. How does the layout spread them across the light and airflow differences?

Compare feeds C and A

The researcher wants C minus A in mean daily gain, in grams per bird per day, adjusted for pen position.

avg_comparisons(
  latin_fit,
  variables = list(feed = c("A", "C")),
  df = df.residual(latin_fit)
)

 Estimate Std. Error    t Pr(>|t|)    S 2.5 % 97.5 % Df
     2.97      0.456 6.52   <0.001 10.7  1.86   4.09  6

Term: feed
Type: response
Comparison: C - A

The model assumes that feed, row and column differences add together. The t interval uses the model’s six remaining error degrees of freedom.

What the square cannot separate

Suppose feed C works especially well only in Row 1. We have one pen at each row-column position.

Can we separately estimate every feed-by-row change, every feed-by-column change, and ordinary pen variation?

Explain why adding more model terms cannot create missing replication.

Unequal sizes are a property of a design

A design is unbalanced when relevant groups do not contain equal numbers of units or measurements.

This can be planned, for example because one treatment is costly, or arise because observations are lost.

Count the pens researchers assigned diets to, as well as the measurements they recorded. Losing animal measurements within pens and losing whole pens are different changes.

Recall the split-plot study

Recall all four rooms

Sixteen pens, eight per diet, 32 section means.

A measuring-equipment fault

A fault leads to the loss of both section averages from pen R1-P01 and the Left section of R3-P02.

Assume this fault is unrelated to animal growth.

u <- read_csv("unequal-split-plot.csv", show_col_types = FALSE) |>
  mutate(
    room = factor(room, levels = paste0("R", 1:4)),
    pen = factor(pen),
    diet = factor(diet, levels = c("Control", "High fibre")),
    bedding = factor(bedding, levels = c("Standard", "Deep"))
  )
u |>
  count(diet, bedding)

# A tibble: 4 × 3
  diet       bedding      n
  <fct>      <fct>    <int>
1 Control    Standard     8
2 Control    Deep         7
3 High fibre Standard     7
4 High fibre Deep         7

Count what remains

There are now fifteen pens with data. Losing a section does not remove its partner’s diet assignment.

remaining_pens <- u |>
  distinct(room, pen, diet)
remaining_pens |>
  count(room, diet)

# A tibble: 8 × 3
  room  diet           n
  <fct> <fct>      <int>
1 R1    Control        2
2 R1    High fibre     1
3 R2    Control        2
4 R2    High fibre     2
5 R3    Control        2
6 R3    High fibre     2
7 R4    Control        2
8 R4    High fibre     2

If the heaviest or sickest animals were more likely to be missing, why would the same count table be insufficient?

Fit the recorded unequal data

Use the same design argument: diet between pens, bedding within pens, and known rooms.

unequal_fit <- lme4::lmer(
  gain_g_day ~ room + diet * bedding + (1 | pen),
  data = u,
  REML = TRUE
)
check_convergence(unequal_fit)

[1] TRUE
attr(,"gradient")
[1] 1.34e-08

check_singularity(unequal_fit)

[1] FALSE

Unequal counts alone do not require deleting good observations to make groups equal.

Raw means count observed sections

These averages give each recorded section equal weight. Rooms with more surviving sections contribute more.

u |>
  group_by(diet, bedding) |>
  summarise(gain_g_day = mean(gain_g_day), .groups = "drop") |>
  arrange(bedding, diet)

# A tibble: 4 × 3
  diet       bedding  gain_g_day
  <fct>      <fct>         <dbl>
1 Control    Standard       747.
2 High fibre Standard       753.
3 Control    Deep           764.
4 High fibre Deep           777.

The farm manager instead wants each of the four rooms to count equally. We need to state that choice before comparing adjusted means.

Estimate means for the same rooms

Use newdata = "balanced" to give each room equal weight within every diet and bedding combination, even after measurements are lost.

avg_predictions(
  unequal_fit,
  newdata = "balanced",
  by = c("diet", "bedding"),
  re.form = NA
)

       diet  bedding Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
 Control    Standard      747       21.2 35.3   <0.001 903.5   706    789
 Control    Deep          762       21.3 35.8   <0.001 929.5   720    804
 High fibre Standard      747       23.0 32.5   <0.001 765.0   702    792
 High fibre Deep          771       23.0 33.5   <0.001 814.7   726    816

Type: response

Each room receives one quarter of the weight. re.form = NA leaves out differences fitted for particular pens; here this also gives the mean over the assumed zero-mean pen differences.

Compare the same target after losses

For each bedding, compare the two diets with equal room weights.

avg_comparisons(
  unequal_fit,
  variables = "diet",
  newdata = "balanced",
  by = "bedding",
  re.form = NA
)

  bedding Estimate Std. Error       z Pr(>|z|)   S 2.5 % 97.5 %
 Standard   -0.225       31.3 -0.0072    0.994 0.0 -61.5   61.1
 Deep        8.928       31.4  0.2848    0.776 0.4 -52.5   70.4

Term: diet
Type: response
Comparison: High fibre - Control

What assumption makes this comparison useful? What information about the fault would you request before using it?

Recall the forty-calf example

These are simulated mean paths, not observed weights or fitted calf results.

The calf random effects follow a normal model

For every calf, we draw a starting-weight random effect from a normal distribution with mean zero and SD 4 kg. We draw a weekly-gain random effect from another normal distribution with mean zero and SD 0.9 kg/week.

We draw each calf’s pair without using another calf’s draw. The pairs are independent and identically distributed, or i.i.d., across calves.

The same pair of distributions applies to every calf. Starting-weight and weekly-gain random effects have different units and spreads. Gaussian means normal.

Keep the random effects in their original units

We want to illustrate the shape of the random-effect distributions. A QQ plot compares their ordered values with a normal distribution; rescaling is unnecessary.

avg_predictions() estimates expected outcomes and avg_comparisons() estimates contrasts. Those answer different questions from checking a distribution’s shape.

We keep starting-weight random effects in kg and weekly-gain random effects in kg/week.

Prepare the calf random effects for plotting

Stack the two columns from the earlier calf illustration. Each panel will show one kind of random effect, in its original units.

random_draws <- calf_effects |>
  transmute(
    calf,
    "Starting weight (kg)" = start_difference,
    "Weekly gain (kg/week)" = gain_difference
  ) |>
  pivot_longer(
    cols = -calf,
    cols_vary = "slowest",
    names_to = "effect",
    values_to = "difference"
  )
random_draws |>
  slice(c(1, 2, 41, 42))

# A tibble: 4 × 3
  calf  effect                difference
  <fct> <chr>                      <dbl>
1 C01   Starting weight (kg)     -13.8  
2 C02   Starting weight (kg)       0.680
3 C01   Weekly gain (kg/week)      1.33 
4 C02   Weekly gain (kg/week)     -1.79 

Build a normal QQ plot

A quantile–quantile plot, or QQ plot, places ordered values against values expected from a reference distribution. Here the reference is normal.

p_random_qq <- ggplot(random_draws, aes(sample = difference)) +
  stat_qq() +
  stat_qq_line(linetype = "dashed") +
  facet_wrap(~effect, nrow = 1, scales = "free_y") +
  labs(
    x = "Quantiles of a standard normal distribution",
    y = "Simulated calf random effect"
  ) +
  theme_bw(base_size = 14)

Each panel keeps its own vertical units. The dashed line provides a normal reference matched to that panel’s values.

Normal QQ plots of the calf random effects

Each panel contains forty simulated random effects. Their original units appear in the panel labels.

p_random_qq

Points near a straight line are consistent with a normal shape. We generated these values from normal distributions; this illustrates the assumption rather than checking real calves.

A QQ plot does not establish independence

Calves sharing a feeder or a disease outbreak might have related random effects. The shape of this QQ plot cannot tell us whether they shared those conditions.

Recorded weights from one calf share its random intercept and slope. We are not assuming that all its repeated weights are separate i.i.d. draws.

After centring time, intercept and slope random effects within a calf can be related while different calves’ pairs still follow the same Gaussian model. Fitted, partly pooled random effects should not be treated as the raw simulated draws shown here.

Return to the split-plot study

Return to all four rooms

Sixteen pens, eight per diet, 32 section means.

Project design and data audit

For your own project, draw where treatment was assigned or observations were sampled.

Mark shared groups, repeated measurements, unequal counts and missing design information.

Record one data source, one unanswered question and who will resolve it before 4 November.

Individual exit

A proposal gives one diet to each pen, and two beddings to sections inside each pen. One section measurement is lost.

Explain why this is both split plot and unbalanced, and why adding more animals would not add independent diet assignments.

After class: a binary example and future power planning with simglm.

Software references

lme4 · Linear models · performance checks · marginaleffects predictions

All numerical examples are invented teaching data.

Backup slides

The greenhouse example is retained for future Agroecology cohorts. The workshop ends before this section.

Backup: fertilisers in a greenhouse

Sixteen trays occupy four rows and four columns. Light changes along the rows; temperature changes along the columns.

We want to compare four fertilisers, A–D, without putting one only in a bright or warm position.

A Latin square places every fertiliser once in each row and once in each column.

Read the square

Each record is one tray’s position, fertiliser and measured mass. A tray receives just one fertiliser.

plant_square <- read_csv("latin-square.csv", show_col_types = FALSE) |>
  mutate(
    row = factor(row, levels = paste0("Row", 1:4)),
    column = factor(column, levels = paste0("Col", 1:4)),
    fertiliser = factor(fertiliser, levels = c("A", "B", "C", "D"))
  )
plant_square |>
  select(row, column, fertiliser) |>
  pivot_wider(names_from = column, values_from = fertiliser)

# A tibble: 4 × 5
  row   Col1  Col2  Col3  Col4 
  <fct> <fct> <fct> <fct> <fct>
1 Row1  D     A     B     C    
2 Row2  A     D     C     B    
3 Row3  C     B     D     A    
4 Row4  B     C     A     D    

Point to all four positions for fertiliser A. Which two sources of variation does this layout spread it across?

Answer the grower’s question

The grower wants C minus A, in grams of plant mass per tray.

avg_comparisons(
  plant_fit,
  variables = list(fertiliser = c("A", "C")),
  df = df.residual(plant_fit)
)

 Estimate Std. Error    t Pr(>|t|)   S 2.5 % 97.5 % Df
     1.38       0.33 4.18  0.00582 7.4 0.571   2.18  6

Term: fertiliser
Type: response
Comparison: C - A

The model assumes the row, column and fertiliser changes add together. The t interval uses six remaining error degrees of freedom.

What the square cannot separate

Suppose fertiliser C works especially well only in Row 1. We have one observation at each row-column position.

Can we separately estimate every fertiliser-by-row change, every fertiliser-by-column change, and ordinary tray variation?

Explain why adding more model terms cannot create missing replication.