Crossed Random Effects in a Penicillin Assay

Linear Mixed Effects Models: Practical Work, Group 5

UNIVERSITY OF ABOMEY-CALAVI (UAC)
Doctoral School of Agronomic and Water Science
Laboratoire de Biomathématiques et d’Estimations Forestières (LaBEF)
Master of Science in Biostatistics, Linear Mixed Effects Models

Submitted to: Dr. Mintode Nicodème ATCHADE

Associate Professor of Statistics and Econometrics · National University of Sciences, Technologies, Engineering and Mathematics (UNSTIM Abomey)

Abstract

Background. In a classical microbiological assay, the diameter of the zone of inhibition produced by six penicillin samples (A to F) was recorded on 24 agar plates, each plate receiving all six samples. Plates and samples are therefore crossed grouping factors, and measurements sharing a plate or a sample are not independent.

Methods. We fitted a linear mixed model with a fixed intercept and crossed random intercepts for plate and sample (REML), compared it with single-factor models (REML and ML, likelihood-ratio tests, AIC and BIC), contrasted it with a naive linear model and with a model where sample is fixed, and checked the assumptions graphically and formally.

Results. The grand mean was 22.97 mm (95% profile CI 21.27 to 24.68). Variance was dominated by differences between samples (3.73, 78.5%), followed by plates (0.72, 15.1%) and residual noise (0.30, 6.4%). The naive model understated the standard error of the mean by a factor of about 4.8 (0.17 against 0.81 mm). Both random terms were strongly supported (LRT, p < 0.001); the fit was not singular and no serious assumption violation was found.

Conclusion. The assay is precise (small residual and plate variances); nearly all the signal comes from true differences between samples. Ignoring the crossed structure would give overconfident inference.

1 Introduction

1.1 Context and rationale

Penicillin potency is traditionally measured by the diffusion assay: a sample is deposited on an agar plate seeded with a susceptible organism, and the diameter of the clear inhibition zone grows with potency. Here, 24 plates each received the six samples, giving 144 measurements (Davies and Goldsmith, 1972; lme4::Penicillin). Two sources of variation compete with the signal: which sample is tested (the effect of interest) and which plate carries it (a nuisance: agar thickness, inoculum density, incubation conditions).

Because every sample appears on every plate, one plate contributes six correlated measurements and one sample contributes 24 correlated measurements. A linear model with independent errors is therefore inappropriate, and a linear mixed model with crossed random effects is the natural tool (Bates et al., 2015).

1.2 Objectives and hypotheses

1. Describe the design and show that plate and sample are crossed and balanced.
2. Estimate and compare the variance contributed by plates, samples and residual error.
3. Show what is lost when the grouping is ignored, and what changes when sample is treated as fixed.
4. Validate the final model and interpret it for a non-statistician.

Hypothesis	Statement
H1	Samples differ in mean inhibition diameter: \(\sigma^2_{sample} > 0\).
H2	Plates contribute non-negligible assay variability: \(\sigma^2_{plate} > 0\).
H3	The crossed model (plate and sample) fits better than either single-factor model.

1.3 Direct answers to the conceptual questions

Q1. Why are plate and sample crossed rather than nested? Model equation. In a nested design each level of one factor occurs within a single level of the other (for example, each cask belongs to one batch). Here every sample is tested on every plate: levels of sample recur across plates and levels of plate recur across samples, so neither is nested in the other (verified on the data in Section 3.1). With \(i = 1,\dots,24\) plates and \(j = 1,\dots,6\) samples:

\[y_{ij} = \mu + a_i + b_j + \varepsilon_{ij}, \qquad a_i \sim N(0,\sigma^2_{a}),\quad b_j \sim N(0,\sigma^2_{b}),\quad \varepsilon_{ij} \sim N(0,\sigma^2),\]

with \(a_i\) (plate), \(b_j\) (sample) and \(\varepsilon_{ij}\) mutually independent. In lme4 syntax: diameter ~ 1 + (1 | plate) + (1 | sample). Implied covariances: \(\mathrm{Var}(y_{ij}) = \sigma^2_a+\sigma^2_b+\sigma^2\); two observations on the same plate share \(\sigma^2_a\); two observations of the same sample share \(\sigma^2_b\); otherwise they are uncorrelated.

Q4. Which assumptions are required for the random effects? (i) Plate effects are independent draws from \(N(0,\sigma_a^2)\) and sample effects from \(N(0,\sigma_b^2)\) (zero mean, normality, constant variance); (ii) plate effects, sample effects and residuals are mutually independent; (iii) residuals are independent, normal and homoscedastic, \(\varepsilon_{ij}\sim N(0,\sigma^2)\); (iv) effects are additive (no plate by sample interaction), because with one observation per cell an interaction cannot be separated from the residual. They are checked in Section 3.8.

General questions (a) and (b). (a) Fixed or random. The intercept \(\mu\) (overall mean) is fixed. Plate is random: the 24 plates are a sample of all plates that could have been prepared, and we want to quantify their variability, not compare them. Sample is random in the prescribed model because the six preparations are regarded as a sample from a wider population and the target is the between-sample variance; the fixed alternative is examined in Section 3.7. (b) Variance components. \(\sigma_a^2\) is plate-to-plate variability (agar, inoculum, incubation: assay noise shared by the six measurements of one plate); \(\sigma_b^2\) is true between-sample variability (potency differences, the signal); \(\sigma^2\) is the residual variability not explained by plate or sample (measurement and local plate noise).

Q3 (principle). When is it more reasonable to treat sample as fixed? When inference concerns these particular six samples (for example “is A more potent than F?”, release of defined batches), when the levels are not exchangeable draws from a population, or when there are too few levels to estimate a variance reliably. Random is justified when the six samples stand for a population and one wants \(\sigma^2_{sample}\) and generalisable inference. The data-based answer is in Section 3.7.

2 Methods

2.1 Data and design

The data frame Penicillin (package lme4) has 144 rows and three variables: diameter (response, mm), plate (factor, 24 levels, a to x) and sample (factor, 6 levels, A to F). Each plate by sample cell holds exactly one observation and there are no missing values.

2.2 Statistical models

All mixed models are fitted with lmerTest::lmer() (Satterthwaite degrees of freedom).

Model	Formula	Purpose
Naive	lm(diameter ~ 1)	Ignores the grouping (shows why it fails)
M1	diameter ~ 1 + (1 \| plate)	Plate variance only
M2	diameter ~ 1 + (1 \| sample)	Sample variance only
M3 (final)	diameter ~ 1 + (1 \| plate) + (1 \| sample)	Crossed random intercepts
M4	diameter ~ sample + (1 \| plate)	Sample as a fixed effect

Under M3 the variance of the grand mean is \(\mathrm{Var}(\bar y_{..}) = \sigma_a^2/24 + \sigma_b^2/6 + \sigma^2/144\): with only six samples, the between-sample variance is divided by 6, not by 144. This is why standard errors computed under independence are far too small.

2.3 Estimation, comparison and inference

• REML estimates and compares variance components (models share the same fixed part). ML is used for comparisons involving different fixed structures (M3 against M4) and to illustrate the small-sample bias of ML variances.
• Random-term tests: likelihood-ratio tests with lmerTest::ranova(). A variance is tested on the boundary of the parameter space (null distribution a 50:50 mixture of \(\chi^2_0\) and \(\chi^2_1\)), so the p-values are conservative; conclusions are also supported by AIC, BIC and profile intervals.
• Uncertainty: profile-likelihood intervals, confint(method = "profile").
• Variance partition and ICC: by hand from VarCorr() and with performance::icc(); \(R^2\) with performance::r2().
• Pairwise comparisons (M4): emmeans with Tukey adjustment.
• Diagnostics: residuals against fitted values, normal Q-Q plots of residuals and random effects, performance::check_model(), outliers, group influence (Cook’s distance, influence.ME), singularity, sensitivity refit. Significance level: 5%.

2.4 Where each requirement is answered

Assignment item	Where
Q1, Q3 (principle), Q4 (assumptions), General (a) and (b)	Introduction
Step 1: design, crossed and balanced, xtabs (To do 1); Q1 verified	3.1
Step 2: graphical exploration	3.2
Step 3: naive model and why it fails	3.3
Steps 4 and 5: model sequence, omitting one random term, REML against ML (To do 2 and 3)	3.4
Step 6: summary, VarCorr, ranef, CIs, ICC, R², dotplot (To do 4); Q2	3.5
General (c): what if the grouping is ignored	3.6
To do 5, Q3: sample as a fixed effect	3.7
Step 7: diagnostics, singularity, influence; Q4 verified	3.8
Step 8: plain-language interpretation; General (d): limitations	4

2.5 Software and reproducibility

needed <- c("lme4", "lmerTest", "performance", "emmeans", "broom", "broom.mixed",
            "lattice", "ggplot2", "dplyr", "tidyr", "tibble", "knitr",
            "kableExtra", "DT", "plotly", "patchwork", "see", "qqplotr",
            "influence.ME", "scales")
to_install <- needed[!needed %in% rownames(installed.packages())]
if (length(to_install) > 0) install.packages(to_install)

library(lme4); library(lmerTest)          # lmerTest last, so lmer() returns p-values
library(ggplot2); library(dplyr); library(tidyr); library(tibble)
library(performance); library(emmeans)
library(patchwork)                        # operators | and / to combine ggplots
library(knitr); library(kableExtra)

set.seed(2026)
pal <- c("#1F5F8B", "#E69F00", "#009E73", "#CC79A7", "#D55E00", "#56B4E9")
theme_set(
  theme_minimal(base_size = 13) +
    theme(plot.title = element_text(face = "bold"),
          plot.subtitle = element_text(colour = "grey35"),
          axis.title = element_text(face = "bold"),
          panel.grid.minor = element_blank(),
          legend.position = "bottom")
)

3 Results

Each output below is preceded by a green box that explains how to read it.

3.1 Design verification

Q1 verified on the data. Each of the six samples is tested on all 24 plates, and the plate by sample table contains exactly one observation in each of the 144 cells. The design is therefore fully crossed, complete and balanced; it is not nested.

data("Penicillin", package = "lme4")
penicillin <- as_tibble(Penicillin)
cross_tab  <- xtabs(~ plate + sample, data = penicillin)

Structure and balance

Interpretation. The first table shows that all 144 possible plate by sample cells are observed once (One obs. per cell = TRUE): complete and balanced. The second table shows that each sample appears on 24 plates, the signature of crossing (in a nested design each sample would appear on a single plate).

tibble(
  Observations        = nrow(penicillin),
  Plates              = nlevels(penicillin$plate),
  Samples             = nlevels(penicillin$sample),
  `Possible cells`    = nlevels(penicillin$plate) * nlevels(penicillin$sample),
  `Observed cells`    = sum(cross_tab > 0),
  `One obs. per cell` = all(cross_tab == 1)
) |>
  kable(caption = "Design summary.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Design summary.

Observations	Plates	Samples	Possible cells	Observed cells	One obs. per cell
144	24	6	144	144	TRUE

penicillin |>
  group_by(Sample = sample) |>
  summarise(`Plates on which it is tested` = n_distinct(plate), .groups = "drop") |>
  kable(caption = "Number of plates per sample.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Number of plates per sample.

Sample	Plates on which it is tested
A	24
B	24
C	24
D	24
E	24
F	24

Layout of the data

Interpretation. Every cell of the grid is filled, which confirms the crossing. Colours change mainly from one column to another (between samples) and much less from one row to another (between plates): high values (dark) sit in columns A and C, low values (light) in column F. This anticipates the variance decomposition of Section 3.5.

ggplot(penicillin, aes(x = reorder(sample, -diameter, FUN = mean),
                       y = reorder(plate, diameter, FUN = mean), fill = diameter)) +
  geom_tile(colour = "white") +
  geom_text(aes(label = diameter), size = 3) +
  scale_fill_gradient(low = "#FFF3D6", high = pal[1], name = "Diameter (mm)") +
  labs(x = "Sample", y = "Plate")

Observed diameters (mm) for each plate by sample cell; plates ordered by mean diameter, samples by mean diameter.

Raw data

Interpretation. The 144 observations are listed so that any value cited in this report can be checked. Use the search boxes to filter by plate or sample.

DT::datatable(penicillin, rownames = FALSE, filter = "top",
              options = list(pageLength = 8, dom = "tip"))

3.2 Exploratory analysis

Interpretation. The 144 diameters average 22.97 mm with an SD of 2.03 mm (median 23, IQR 2, range 18 to 27). The mean and median coincide, which suggests a roughly symmetric marginal distribution.

penicillin |>
  summarise(n = n(), Mean = mean(diameter), SD = sd(diameter),
            Median = median(diameter), IQR = IQR(diameter),
            Min = min(diameter), Max = max(diameter)) |>
  kable(digits = 2, caption = "Overall descriptive statistics of the inhibition diameter (mm).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Overall descriptive statistics of the inhibition diameter (mm).

n	Mean	SD	Median	IQR	Min	Max
144	22.97	2.03	23	2	18	27

Interpretation. Sample means range from 19.96 mm (F) to 25.17 mm (A), a spread of about 5 mm, whereas the SD within each sample is only about 1 mm (0.95 to 1.08). Most of the total variability (SD 2.03) is therefore between samples. Samples A and C are the most potent, F the least.

sample_summary <- penicillin |>
  group_by(Sample = sample) |>
  summarise(n = n(), Mean = mean(diameter), SD = sd(diameter),
            Min = min(diameter), Max = max(diameter), .groups = "drop") |>
  arrange(desc(Mean))
kable(sample_summary, digits = 2, caption = "Descriptive statistics by sample, sorted by mean.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Descriptive statistics by sample, sorted by mean.

Sample	n	Mean	SD	Min	Max
A	24	25.17	0.96	24	27
C	24	24.92	1.02	23	26
E	24	22.96	1.04	21	25
D	24	22.88	0.95	21	24
B	24	21.96	1.00	20	24
F	24	19.96	1.08	18	22

By sample

Interpretation. Boxes of A and C sit around 25 mm, those of D and E around 23 mm, B around 22 mm and F around 20 mm. Within a sample the points are tightly grouped (about 1 mm), so the boxes of A and C do not overlap with those of F: samples clearly differ. Hover over a point to read its value.

p_sample <- ggplot(penicillin, aes(x = reorder(sample, -diameter, FUN = mean),
                                   y = diameter, fill = sample)) +
  geom_boxplot(width = 0.6, alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.1, alpha = 0.6, size = 1.6) +
  stat_summary(fun = mean, geom = "point", shape = 23, size = 3, fill = "white") +
  scale_fill_manual(values = pal, guide = "none") +
  labs(title = "Inhibition diameter by sample",
       subtitle = "Diamonds: sample means; points: individual plates",
       x = "Sample (ordered by mean)", y = "Diameter (mm)")
plotly::ggplotly(p_sample)

By plate

Interpretation. Plate means vary between roughly 21.5 and 24.5 mm (plate m is the highest), a much narrower range than the differences between samples. Each box is wide because it contains the six samples of that plate: the width of the boxes reflects the sample effect, while the shifts between boxes reflect the plate effect.

p_plate <- ggplot(penicillin, aes(x = reorder(plate, diameter, FUN = mean),
                                  y = diameter)) +
  geom_boxplot(fill = pal[3], alpha = 0.6, outlier.shape = NA) +
  geom_jitter(width = 0.1, alpha = 0.6, size = 1.3) +
  labs(title = "Inhibition diameter by plate",
       subtitle = "Plates ordered by mean diameter",
       x = "Plate", y = "Diameter (mm)")
plotly::ggplotly(p_plate)

Plate profiles

Interpretation. Each line joins the six measurements of one plate. The lines are close to parallel and follow the same ordering of the samples (high at A and C, low at F): sample effects are consistent across plates, and a plate mostly moves its whole profile up or down. This supports the additive structure of the model (no strong plate by sample interaction).

p_prof <- ggplot(penicillin, aes(x = sample, y = diameter, group = plate, colour = plate)) +
  geom_line(alpha = 0.6) + geom_point(size = 1.4) +
  labs(title = "Plate-specific profiles across samples",
       x = "Sample", y = "Diameter (mm)") +
  theme(legend.position = "none")
plotly::ggplotly(p_prof, tooltip = c("x", "y", "colour"))

Distribution

Interpretation. The marginal distribution is unimodal, fairly symmetric (slightly longer on the left) and has no extreme value, with most diameters between 21 and 26 mm. The mixed model assumes normality of residuals and random effects, not of the raw response; this plot only shows that the response is well behaved.

ggplot(penicillin, aes(x = diameter)) +
  geom_histogram(aes(y = after_stat(density)), bins = 15,
                 fill = pal[1], colour = "white", alpha = 0.85) +
  geom_density(colour = pal[5], linewidth = 1) +
  labs(x = "Diameter (mm)", y = "Density")

Histogram and density of the inhibition diameter.

3.3 Naive model ignoring the grouping

m_lm0       <- lm(diameter ~ 1, data = penicillin)
m_lm_sample <- lm(diameter ~ sample, data = penicillin)
m_lm_both   <- lm(diameter ~ plate + sample, data = penicillin)

Interpretation. If the observations were independent, the residuals of the intercept-only model would look like random noise around zero. Instead they are organised by sample: residuals of A and C are about +2 mm, those of F about -3 mm, those of B about -1 mm. Residuals of the same sample share a common offset, so independence is violated and lm() is inappropriate.

tibble(Sample = penicillin$sample, Residual = resid(m_lm0)) |>
  ggplot(aes(x = Sample, y = Residual, fill = Sample)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.1, alpha = 0.5) +
  scale_fill_manual(values = pal, guide = "none") +
  labs(y = "Residual (mm)")

Residuals of lm(diameter ~ 1), grouped by sample.

Interpretation. The residual SD drops from 2.03 mm (nothing modelled) to about 1.0 mm once sample is included, and to about 0.55 mm once plate is also included. Most of what the naive model calls “noise” is in fact sample and plate structure; the correct standard errors are compared in Section 3.6.

tibble(
  Model = c("lm(diameter ~ 1)", "lm(diameter ~ sample)", "lm(diameter ~ plate + sample)"),
  `Residual SD (mm)` = c(sigma(m_lm0), sigma(m_lm_sample), sigma(m_lm_both)),
  `Residual df` = c(df.residual(m_lm0), df.residual(m_lm_sample), df.residual(m_lm_both))
) |>
  kable(digits = 3, caption = "Residual SD of fixed-effects-only models.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Residual SD of fixed-effects-only models.

Model	Residual SD (mm)	Residual df
lm(diameter ~ 1)	2.031	143
lm(diameter ~ sample)	1.010	138
lm(diameter ~ plate + sample)	0.550	115

3.4 Mixed-model building and comparison

m_plate  <- lmer(diameter ~ 1 + (1 | plate),                data = penicillin)
m_sample <- lmer(diameter ~ 1 + (1 | sample),               data = penicillin)
m_full   <- lmer(diameter ~ 1 + (1 | plate) + (1 | sample), data = penicillin)
m_final  <- m_full

m_plate_ml  <- update(m_plate,  REML = FALSE)
m_sample_ml <- update(m_sample, REML = FALSE)
m_full_ml   <- update(m_full,   REML = FALSE)

Information criteria

Interpretation. The crossed model M3 has by far the lowest AIC and BIC (REML: AIC 338.9 against 441.9 for sample only and 619.3 for plate only; ML gives the same ranking). A difference of more than 100 points is overwhelming evidence that both grouping factors are needed (hypothesis H3). Note that the sample term matters much more than the plate term: dropping sample costs 280 AIC points, dropping plate costs 103.

lab <- c("Plate only (M1)", "Sample only (M2)", "Plate + Sample (M3)")
ic_reml <- AIC(m_plate, m_sample, m_full) |>
  mutate(BIC = BIC(m_plate, m_sample, m_full)$BIC, Model = lab, Fit = "REML")
ic_ml <- AIC(m_plate_ml, m_sample_ml, m_full_ml) |>
  mutate(BIC = BIC(m_plate_ml, m_sample_ml, m_full_ml)$BIC, Model = lab, Fit = "ML")
bind_rows(ic_reml, ic_ml) |>
  select(Fit, Model, Parameters = df, AIC, BIC) |>
  kable(digits = 1, row.names = FALSE, caption = "AIC and BIC of the three random-effects structures.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) |>
  row_spec(c(3, 6), bold = TRUE, background = "#E8F6F1")

AIC and BIC of the three random-effects structures.

Fit	Model	Parameters	AIC	BIC
REML	Plate only (M1)	3	619.3	628.2
REML	Sample only (M2)	3	441.9	450.8
REML	Plate + Sample (M3)	4	338.9	350.7
ML	Plate only (M1)	3	617.6	626.5
ML	Sample only (M2)	3	443.2	452.1
ML	Plate + Sample (M3)	4	340.2	352.1

Likelihood-ratio tests

Interpretation. Each row tests one random term by removing it from M3. Removing plate gives \(\chi^2 = 105.1\) and removing sample gives \(\chi^2 = 282.4\) (1 df each), both with p < 0.001. Because a variance is tested at the edge of its parameter space, these p-values are conservative, which only strengthens the conclusion: H1 and H2 are supported, and both variance components are non-zero.

ranova(m_full) |>
  as.data.frame() |>
  rownames_to_column("Term dropped") |>
  mutate(`Pr(>Chisq)` = ifelse(is.na(`Pr(>Chisq)`), "",
                               format.pval(`Pr(>Chisq)`, digits = 2, eps = 0.001))) |>
  kable(digits = 2, caption = "Likelihood-ratio test of each random term (REML).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Likelihood-ratio test of each random term (REML).

Term dropped	npar	logLik	AIC	LRT	Df	Pr(>Chisq)
<none>	4	-165.4	338.9
(1 &#124; plate)	3	-218.0	441.9	105.1	1	<0.001
(1 &#124; sample)	3	-306.6	619.3	282.4	1	<0.001

REML against ML

Interpretation. Plate and residual variances are identical under both methods, but ML gives a sample variance of 3.14 against 3.73 under REML (a ratio of 0.84, about 16% too small). ML ignores the degree of freedom used to estimate the mean, which matters when the grouping factor has only six levels. REML is therefore preferred for variance components here, and ML is kept for comparing fixed structures.

left_join(
  as.data.frame(VarCorr(m_full))    |> select(Component = grp, REML = vcov),
  as.data.frame(VarCorr(m_full_ml)) |> select(Component = grp, ML = vcov),
  by = "Component"
) |>
  mutate(`ML / REML` = ML / REML) |>
  kable(digits = 3, caption = "Variance components under REML and ML.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Variance components under REML and ML.

Component	REML	ML	ML / REML
plate	0.717	0.715	0.997
sample	3.731	3.135	0.840
Residual	0.302	0.302	1.000

3.5 Final model inspection

ci_prof <- confint(m_final, method = "profile", oldNames = FALSE)

Summary

Interpretation. The model converged without warning on 144 observations, 24 plates and 6 samples. The fixed intercept (grand mean) is 22.97 mm with a standard error of 0.81 mm and only 5.5 degrees of freedom (an effective count driven by the six samples); it is highly significant (t = 28.4, p < 0.001), which only says that the mean diameter is not zero. The random-effects block gives the SDs: plate 0.85 mm, sample 1.93 mm, residual 0.55 mm. Scaled residuals lie between -2.08 and 2.98, a reasonable range for 144 points.

summary(m_final)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: diameter ~ 1 + (1 | plate) + (1 | sample)
##    Data: penicillin
## 
## REML criterion at convergence: 330.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0792 -0.6714  0.0629  0.5838  2.9796 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plate    (Intercept) 0.717    0.847   
##  sample   (Intercept) 3.731    1.932   
##  Residual             0.302    0.550   
## Number of obs: 144, groups:  plate, 24; sample, 6
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   22.972      0.809  5.487    28.4  3.6e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Variance components (Q2)

Q2. Which variance is larger, between-plate or between-sample? What does it mean for the assay? The between-sample variance (3.73, SD 1.93 mm, 78.5% of the total) is about 5.2 times larger than the between-plate variance (0.72, SD 0.85 mm, 15.1%); the residual variance is small (0.30, 6.4%). Meaning: the assay does what it is designed to do, since the dominant source of variation is the true difference between samples (the signal). Plate-to-plate variation is secondary but real (the LRT is highly significant), which justifies testing all samples on every plate so that plate effects cancel in comparisons between samples. The small residual variance indicates good repeatability. The sample SD is, however, imprecise (95% CI 1.10 to 3.56 mm) because there are only six samples.

Interpretation. The table converts the three variances into shares of the total: 78.5% sample, 15.1% plate, 6.4% residual. The bar chart shows the same decomposition visually: the sample bar dominates, the plate bar is about five times shorter, and the residual is small.

vc_tab <- as.data.frame(VarCorr(m_final)) |>
  transmute(Component = recode(grp, plate = "Plate", sample = "Sample"),
            Variance = vcov, SD = sdcor) |>
  mutate(`% of total` = 100 * Variance / sum(Variance))
kable(vc_tab, digits = 3, caption = "Variance components of the final model (REML).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Variance components of the final model (REML).

Component	Variance	SD	% of total
Plate	0.717	0.847	15.091
Sample	3.731	1.932	78.543
Residual	0.302	0.550	6.366

ggplot(vc_tab, aes(x = reorder(Component, `% of total`), y = `% of total`, fill = Component)) +
  geom_col(width = 0.6, show.legend = FALSE) +
  geom_text(aes(label = sprintf("%.1f%%", `% of total`)), hjust = -0.15, fontface = "bold") +
  coord_flip(ylim = c(0, 95)) +
  scale_fill_manual(values = c(Plate = pal[2], Sample = pal[1], Residual = pal[4])) +
  labs(x = NULL, y = "Share of total variance (%)")

Share of the total variance attributable to each component (REML).

Estimates and CIs

Interpretation. With 95% confidence, the grand mean lies between 21.27 and 24.68 mm. The plate SD (0.85 mm) is estimated fairly precisely (0.63 to 1.18), the residual SD (0.55 mm) very precisely (0.49 to 0.63), but the sample SD (1.93 mm) is very uncertain (1.10 to 3.56) since it rests on six levels. None of the SD intervals includes zero, which agrees with the LRTs.

tibble(
  Parameter = c("SD plate", "SD sample", "SD residual", "Intercept (grand mean, mm)"),
  Estimate  = c(sqrt(vc_tab$Variance[1]), sqrt(vc_tab$Variance[2]),
                sqrt(vc_tab$Variance[3]), fixef(m_final)[["(Intercept)"]]),
  `2.5 %`   = ci_prof[, 1],
  `97.5 %`  = ci_prof[, 2]
) |>
  kable(digits = 3, caption = "Estimates with 95% profile-likelihood confidence intervals.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Estimates with 95% profile-likelihood confidence intervals.

Parameter	Estimate	2.5 %	97.5 %
SD plate	0.847	0.634	1.182
SD sample	1.932	1.096	3.556
SD residual	0.550	0.486	0.629
Intercept (grand mean, mm)	22.972	21.267	24.678

Interpretation. The t-test of the intercept (t = 28.4 on about 5.5 df, p < 0.001) confirms that the mean diameter differs from zero. This test is not of scientific interest; the useful quantities are the variance components above.

coef(summary(m_final)) |>
  as.data.frame() |>
  mutate(`Pr(>|t|)` = format.pval(`Pr(>|t|)`, digits = 2, eps = 0.001)) |>
  kable(digits = 3, caption = "Test of the fixed intercept (Satterthwaite df).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Test of the fixed intercept (Satterthwaite df).

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	22.97	0.809	5.487	28.41	<0.001

ICC and R2

Interpretation. The intraclass correlation of each factor is its share of the total variance: sample 0.785 and plate 0.151. Two measurements of the same sample on different plates are strongly correlated (ICC 0.785), which is why independence fails. The two factors together explain 93.6% of the variance (adjusted ICC 0.936), leaving only 6.4% as residual noise.

vc_tab |>
  transmute(Component, ICC = Variance / sum(Variance)) |>
  kable(digits = 3, caption = "Intraclass correlation by hand: variance share per component.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Intraclass correlation by hand: variance share per component.

Component	ICC
Plate	0.151
Sample	0.785
Residual	0.064

as.data.frame(performance::icc(m_final)) |>
  kable(digits = 3, caption = "performance::icc(): adjusted and unadjusted ICC (both factors together).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

performance::icc(): adjusted and unadjusted ICC (both factors together).

ICC_adjusted	ICC_unadjusted	optional
0.936	0.936	FALSE

as.data.frame(performance::icc(m_final, by_group = TRUE)) |>
  kable(digits = 3, caption = "performance::icc(by_group = TRUE): ICC of each factor.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

performance::icc(by_group = TRUE): ICC of each factor.

Group	ICC
plate	0.151
sample	0.785

Interpretation. Marginal R² is 0 because the model has no fixed predictor, and conditional R² is 0.936: the random effects account for 93.6% of the variance. The residual SD (sigma) is 0.55 mm and the RMSE 0.495 mm, which is small compared with a total SD of 2.03 mm.

performance::r2(m_final)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.936
##      Marginal R2: 0.000

performance::model_performance(m_final) |>
  kable(digits = 3, caption = "Global performance indices.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Global performance indices.

AIC	AICc	BIC	R2_conditional	R2_marginal	ICC	RMSE	Sigma
338.9	339.1	350.7	0.936	0	0.936	0.495	0.55

Random effects: plates

Interpretation. Each point is the predicted deviation of a plate from the grand mean, with its 95% interval. Plates range from about +1.4 mm (m) to about -1.4 mm (g and s). Several plates (m, t, o, k high; x, s, g low) have intervals that exclude zero, so some plates really give larger or smaller zones than average, which is the plate variance made visible. Many other plates are close to zero.

re_dot <- lattice::dotplot(ranef(m_final, condVar = TRUE))
print(re_dot$plate)

lattice::dotplot of plate random intercepts with conditional intervals.

Random effects: samples

Interpretation. Sample effects are much larger than plate effects: A (+2.19 mm) and C (+1.94 mm) are well above average, F (-3.00 mm) and B (-1.01 mm) well below, while D and E are close to the average. The intervals for A, B, C and F exclude zero. These are the differences the assay is designed to detect.

print(re_dot$sample)

lattice::dotplot of sample random intercepts with conditional intervals.

Interactive caterpillar plot

Interpretation. This is the same information as the two dot plots, on a common scale and interactive (hover to read each estimate and interval). The scale makes the contrast obvious: the sample effects spread over about 5 mm, the plate effects over about 3 mm, although the former are only six points.

re_tab <- broom.mixed::tidy(m_final, effects = "ran_vals", conf.int = TRUE)
p_re <- ggplot(re_tab, aes(x = estimate, y = reorder(level, estimate), colour = group)) +
  geom_vline(xintercept = 0, linetype = 2, colour = "grey50") +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size = 0.35) +
  facet_wrap(~ group, scales = "free_y") +
  scale_colour_manual(values = c(plate = pal[2], sample = pal[1]), guide = "none") +
  labs(title = "Predicted random intercepts with 95% intervals",
       x = "Deviation from the grand mean (mm)", y = NULL)
plotly::ggplotly(p_re)

Shrinkage

Interpretation. Shrinkage pulls predicted effects slightly towards zero. Here it is negligible (at most 0.010 mm, for sample F), because each sample is measured 24 times, so its raw mean is already very reliable. The raw and the model-based sample deviations are practically identical.

grand_mean  <- mean(penicillin$diameter)
blup_sample <- re_tab |> filter(group == "sample") |> select(Sample = level, BLUP = estimate)
penicillin |>
  group_by(Sample = as.character(sample)) |>
  summarise(`Raw deviation` = mean(diameter) - grand_mean, .groups = "drop") |>
  left_join(blup_sample, by = "Sample") |>
  mutate(Shrinkage = `Raw deviation` - BLUP) |>
  kable(digits = 3, caption = "Sample effects: raw deviation from the grand mean against BLUP (mm).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Sample effects: raw deviation from the grand mean against BLUP (mm).

Sample	Raw deviation	BLUP	Shrinkage
A	2.194	2.187	0.007
B	-1.014	-1.010	-0.003
C	1.944	1.938	0.007
D	-0.097	-0.097	0.000
E	-0.014	-0.014	0.000
F	-3.014	-3.004	-0.010

3.6 What if the grouping is ignored?

General (c). How would the conclusions change if the grouping structure were ignored? The naive model (i) treats 144 correlated measurements as independent and reports a standard error of 0.17 mm instead of 0.81 mm (about 4.8 times too small), hence a confidence interval of about 22.6 to 23.3 mm instead of 21.3 to 24.7 mm; (ii) behaves as if it had 144 independent observations whereas the effective sample size for the mean is only about 7; (iii) assigns 93.6% of the variance to “random noise” (residual SD 2.03 mm instead of 0.55 mm), hiding the plate and sample structure. Tests and intervals would be overconfident, with an inflated type I error.

Interpretation. Both models give the same estimate (22.97 mm), but the uncertainty differs by a factor of 4.8. The mixed-model standard error (0.809 mm) is recovered exactly from \(\sigma_a^2/24+\sigma_b^2/6+\sigma^2/144\) (variance 0.654), and the ratio of total variance to variance of the mean gives an effective sample size of about 7.3 independent observations, out of 144 measurements.

vc <- as.data.frame(VarCorr(m_final))
v_plate  <- vc$vcov[vc$grp == "plate"]
v_sample <- vc$vcov[vc$grp == "sample"]
v_resid  <- vc$vcov[vc$grp == "Residual"]

var_mean <- v_plate / 24 + v_sample / 6 + v_resid / 144    # variance of the grand mean
n_eff    <- (v_plate + v_sample + v_resid) / var_mean      # effective sample size

naive_ci <- confint(m_lm0)["(Intercept)", ]
mixed_ci <- ci_prof["(Intercept)", ]
tibble(
  Model = c("Naive lm (independence)", "Crossed mixed model"),
  Estimate = c(coef(m_lm0)[[1]], fixef(m_final)[[1]]),
  `Std. error` = c(coef(summary(m_lm0))[1, 2], coef(summary(m_final))[1, 2]),
  `CI lower` = c(naive_ci[[1]], mixed_ci[[1]]),
  `CI upper` = c(naive_ci[[2]], mixed_ci[[2]])
) |>
  kable(digits = 3, caption = "Grand mean (mm): naive against mixed-model uncertainty.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Grand mean (mm): naive against mixed-model uncertainty.

Model	Estimate	Std. error	CI lower	CI upper
Naive lm (independence)	22.97	0.169	22.64	23.31
Crossed mixed model	22.97	0.809	21.27	24.68

tibble(`Variance of the mean` = var_mean, `Effective sample size` = n_eff) |>
  kable(digits = 3, caption = "Variance of the grand mean and effective sample size under M3.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Variance of the grand mean and effective sample size under M3.

Variance of the mean	Effective sample size
0.654	7.266

3.7 Sample as a fixed effect

Q3. When is it more reasonable to treat sample as fixed? What changes in the inference? Sample should be fixed when interest lies in these six preparations. Then (M4): the F-test for sample is overwhelming (F = 297 on 5 and 115 df); five contrasts are estimated against sample A (B, D, E and F are 2.2 to 5.2 mm smaller; C is not different, p = 0.11); and Tukey comparisons show that all pairs differ except A and C, and D and E. The intercept becomes the mean of sample A (25.17 mm, SE 0.20) instead of a population mean (22.97 mm, SE 0.81), and conclusions apply to these six samples only. Treating sample as random instead yields a variance component and a population-level mean with an honest, much larger, standard error. Point estimates per sample are nearly identical in both approaches; what changes is the question answered and the uncertainty of the overall mean.

m_sfix    <- lmer(diameter ~ sample + (1 | plate), data = penicillin)   # REML
m_sfix_ml <- update(m_sfix, REML = FALSE)
emm       <- emmeans(m_sfix, ~ sample)

F-test and means

Interpretation. The F-test rejects equality of the six sample means (F = 297.1, 5 and 115 df, p < 0.001). The estimated marginal means (with SE 0.21 mm each) reproduce the descriptive means: A 25.17, C 24.92, E 22.96, D 22.87, B 21.96, F 19.96 mm.

anova(m_sfix) |>
  as.data.frame() |>
  mutate(`Pr(>F)` = format.pval(`Pr(>F)`, digits = 2, eps = 0.001)) |>
  kable(digits = 2, caption = "F-test for sample (Satterthwaite df).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

F-test for sample (Satterthwaite df).

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
sample	449.2	89.84	5	115	297.1	<0.001

as.data.frame(emm) |>
  kable(digits = 2, caption = "Estimated marginal means per sample (mm).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Estimated marginal means per sample (mm).

sample	emmean	SE	df	lower.CL	upper.CL
A	25.17	0.21	39.73	24.75	25.58
B	21.96	0.21	39.73	21.54	22.37
C	24.92	0.21	39.73	24.50	25.33
D	22.87	0.21	39.73	22.46	23.29
E	22.96	0.21	39.73	22.54	23.37
F	19.96	0.21	39.73	19.54	20.37

Pairwise comparisons

Interpretation. Of the 15 Tukey-adjusted comparisons, 13 are significant (p < 0.001). The only exceptions are A against C (difference 0.25 mm, p = 0.62) and D against E (difference 0.08 mm, p = 0.995). The most potent group (A, C), the intermediate pair (D, E), B, and F therefore form four distinguishable levels. The largest gap is A against F (5.2 mm).

as.data.frame(pairs(emm, adjust = "tukey")) |>
  mutate(Significant = ifelse(p.value < 0.05, "yes", "no"),
         p.value = format.pval(p.value, digits = 2, eps = 0.001)) |>
  kable(digits = 3, caption = "Tukey-adjusted pairwise differences between samples (mm).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Tukey-adjusted pairwise differences between samples (mm).

contrast	estimate	SE	df	t.ratio	p.value	Significant
A - B	3.208	0.159	115	20.210	<0.001	yes
A - C	0.250	0.159	115	1.575	0.62	no
A - D	2.292	0.159	115	14.436	<0.001	yes
A - E	2.208	0.159	115	13.911	<0.001	yes
A - F	5.208	0.159	115	32.809	<0.001	yes
B - C	-2.958	0.159	115	-18.635	<0.001	yes
B - D	-0.917	0.159	115	-5.774	<0.001	yes
B - E	-1.000	0.159	115	-6.299	<0.001	yes
B - F	2.000	0.159	115	12.598	<0.001	yes
C - D	2.042	0.159	115	12.861	<0.001	yes
C - E	1.958	0.159	115	12.336	<0.001	yes
C - F	4.958	0.159	115	31.234	<0.001	yes
D - E	-0.083	0.159	115	-0.525	1.00	no
D - F	2.917	0.159	115	18.373	<0.001	yes
E - F	3.000	0.159	115	18.898	<0.001	yes

Plot of the comparisons

Interpretation. Each blue bar is the confidence interval of a sample mean and the red arrows are the Tukey comparison arrows: two means differ significantly when their arrows do not overlap. Arrows overlap only for A and C and for D and E, which matches the table above.

plot(emm, comparisons = TRUE) +
  labs(x = "Estimated mean diameter (mm)", y = "Sample")

Estimated mean diameter per sample with Tukey comparison arrows.

Fixed against random

Interpretation. The plate SD (0.85 against 0.83) and residual SD (0.55 against 0.54) barely change between the two formulations. What changes is the intercept (population mean against mean of sample A), its standard error (0.81 against 0.20) and what is estimated for samples (one variance against five contrasts). M4 has the lower AIC (311.7 against 340.2), but this does not settle the choice: the two models answer different questions (these six samples against a population of samples), so the decision must follow the scientific target.

tibble(
  Aspect = c("Target of inference", "Intercept", "SE of the intercept",
             "Sample information", "SD plate", "SD residual", "AIC (ML)"),
  `Sample random (M3)` = c("Population of samples; variance component",
                            "22.97 = grand mean", "0.81 (df about 5.5)",
                            "One SD: 1.93 (CI 1.10 to 3.56)",
                            "0.85", "0.55", format(round(AIC(m_full_ml), 1), nsmall = 1)),
  `Sample fixed (M4)` = c("These six samples; contrasts",
                           "25.17 = mean of sample A", "0.20",
                           "Five contrasts against A (range 5.2 mm)",
                           "0.83", "0.54", format(round(AIC(m_sfix_ml), 1), nsmall = 1))
) |>
  kable(caption = "What changes in the inference when sample is fixed.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = TRUE)

What changes in the inference when sample is fixed.

Aspect	Sample random (M3)	Sample fixed (M4)
Target of inference	Population of samples; variance component	These six samples; contrasts
Intercept	22.97 = grand mean	25.17 = mean of sample A
SE of the intercept	0.81 (df about 5.5)	0.20
Sample information	One SD: 1.93 (CI 1.10 to 3.56)	Five contrasts against A (range 5.2 mm)
SD plate	0.85	0.83
SD residual	0.55	0.54
AIC (ML)	340.2	311.7

3.8 Model diagnostics

Q4 verified on the data. Residuals show no curvature or funnel shape and no evidence against normality (p = 0.30) or constant variance (p = 0.97); the random-effects Q-Q plots are compatible with normality as far as 24 plates and 6 samples allow; the fit is not singular; no plate is influential. Sample F is flagged by Cook’s distance, but only because it is the most extreme real sample (see the influence tab). The independence and additivity assumptions follow from the design (each cell measured once) and from the near-parallel profiles of Section 3.2.

diag_df <- tibble(
  Fitted = fitted(m_final),
  Residual = resid(m_final),
  Scaled = resid(m_final, scaled = TRUE),
  Plate = penicillin$plate,
  Sample = penicillin$sample,
  Obs = seq_len(nrow(penicillin))
)

Residual plots

Interpretation. Residuals against fitted values: the smoother is flat around zero and the vertical spread is constant, so linearity and homoscedasticity hold. The diagonal streaks come from the measurements being rounded to whole millimetres, not from a model defect. Q-Q plot of residuals: points follow the line, with only slightly light lower tail and one high point (observation 137). Q-Q plots of random effects: plate effects are close to the line, with small departures at both ends; the six sample effects follow a roughly linear pattern, but six points cannot really test normality.

p1 <- ggplot(diag_df, aes(Fitted, Residual)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "loess", se = FALSE, colour = pal[5], linewidth = 0.9) +
  labs(title = "Residuals against fitted", x = "Fitted (mm)", y = "Conditional residual")
p2 <- ggplot(diag_df, aes(sample = Scaled)) +
  stat_qq(alpha = 0.7) + stat_qq_line(colour = pal[5]) +
  labs(title = "Normal Q-Q of residuals", x = "Theoretical", y = "Scaled residual")
p3 <- ggplot(re_tab, aes(sample = estimate)) +
  stat_qq(alpha = 0.8) + stat_qq_line(colour = pal[5]) +
  facet_wrap(~ group, scales = "free") +
  labs(title = "Normal Q-Q of random effects",
       subtitle = "Only 6 sample levels: interpret with caution",
       x = "Theoretical", y = "Estimated effect")
(p1 | p2) / p3

Residual diagnostics of the final model.

performance::check_model

Interpretation. The posterior predictive check shows that data simulated from the model (blue) reproduce the observed density (green), with a small discrepancy on the right shoulder. The linearity and homogeneity panels have flat reference lines. The normality panel stays within its confidence band, and the random-effects panels (plate and sample) lie along their reference lines within the bands. No panel signals a problem.

check_model(m_final, check = c("pp_check", "linearity", "homogeneity", "qq", "reqq"))

performance::check_model() panels for the final model.

Formal checks

Interpretation. The normality test of residuals (p = 0.300) and the test of constant variance (p = 0.969) are both far above 0.05, so neither assumption is rejected. check_singularity() and isSingular() return FALSE: all variance components are clearly positive and the fit is not on the boundary. Formal tests complement, but do not replace, the plots above.

check_normality(m_final)

## OK: residuals appear as normally distributed (p = 0.300).

check_heteroscedasticity(m_final)

## OK: Error variance appears to be homoscedastic (p = 0.969).

check_singularity(m_final)

## [1] FALSE

isSingular(m_final, tol = 1e-5)

## [1] FALSE

Outliers

Interpretation. Only two observations have a scaled residual beyond 2 in absolute value: observation 137 (plate w, sample E, +2.98, about three SDs) and observation 14 (plate c, sample B, -2.08). With 144 observations, about seven values beyond |2| would be expected under normality, so this is, if anything, fewer than expected. Observation 137 is the only mildly extreme point; it is kept, and its effect is checked in the sensitivity tab.

diag_df |>
  filter(abs(Scaled) > 2) |>
  arrange(desc(abs(Scaled))) |>
  select(Obs, Plate, Sample, Fitted, Residual, Scaled) |>
  kable(digits = 2, caption = "Observations with |scaled residual| above 2.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Observations with |scaled residual| above 2.

Obs	Plate	Sample	Fitted	Residual	Scaled
137	w	E	22.36	1.64	2.98
14	c	B	22.14	-1.14	-2.08

Influential groups

Interpretation. Cook’s distance measures how much the estimates change when a whole group is removed. Plates: the largest value (plate m) is about 0.007, far below the usual 4/24 = 0.17 threshold, so no plate is influential. Samples: sample F reaches about 0.87, above the crude 4/6 = 0.67 threshold, and is flagged in red. This is expected, since F is the most extreme sample (3 mm below the mean) and with only six samples each one carries much weight on the sample variance. F is a genuine sample, not a measurement error, so it is kept; the consequence is that the sample variance estimate depends strongly on F.

library(influence.ME)
m_lme4      <- lme4::lmer(diameter ~ 1 + (1 | plate) + (1 | sample), data = penicillin)
infl_plate  <- influence(m_lme4, group = "plate")
infl_sample <- influence(m_lme4, group = "sample")
par(mfrow = c(1, 2))
plot(infl_plate,  which = "cook", cutoff = 4 / 24, sort = TRUE,
     xlab = "Cook's distance", ylab = "Plate", main = "Influence of plates")

Cook’s distance of each plate (left) and each sample (right).

plot(infl_sample, which = "cook", cutoff = 4 / 6, sort = TRUE,
     xlab = "Cook's distance", ylab = "Sample", main = "Influence of samples")

Cook’s distance of each plate (left) and each sample (right).

par(mfrow = c(1, 1))

Sensitivity refit

Interpretation. Removing the most extreme observation (plate w, sample E) leaves the conclusions unchanged: the sample variance is 3.731 against 3.734, the plate variance moves from 0.717 to 0.746, the residual variance from 0.302 to 0.274, and the grand mean from 22.972 to 22.958 mm. The model is therefore robust to this observation, which is retained.

m_sens <- update(m_final, data = filter(penicillin, !(plate == "w" & sample == "E")))
bind_rows(
  Full    = as.data.frame(VarCorr(m_final)),
  Without = as.data.frame(VarCorr(m_sens)),
  .id = "Dataset"
) |>
  select(Dataset, Component = grp, vcov) |>
  pivot_wider(names_from = Dataset, values_from = vcov) |>
  bind_rows(tibble(Component = "Intercept", Full = fixef(m_final)[[1]], Without = fixef(m_sens)[[1]])) |>
  kable(digits = 3, caption = "Estimates with and without observation 137 (plate w, sample E).") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Estimates with and without observation 137 (plate w, sample E).

Component	Full	Without
plate	0.717	0.746
sample	3.731	3.734
Residual	0.302	0.274
Intercept	22.972	22.958

Observed against fitted

Interpretation. Conditional predictions (which include the plate and sample effects) fall close to the diagonal for all six samples, with a typical error of about 0.5 mm (RMSE 0.495). The samples form six distinct clouds along the diagonal, a visual reminder that the sample effect explains most of the variation. This is consistent with the conditional R² of 0.936.

ggplot(penicillin |> mutate(Fitted = predict(m_final)),
       aes(x = Fitted, y = diameter, colour = sample)) +
  geom_abline(intercept = 0, slope = 1, linetype = 2, colour = "grey40") +
  geom_point(size = 2, alpha = 0.8) +
  scale_colour_manual(values = pal, name = "Sample") +
  labs(x = "Fitted (mm)", y = "Observed (mm)")

Observed against conditional fitted diameters; dashed line: perfect agreement.

4 Discussion

4.1 Plain-language interpretation

Imagine the six penicillin samples being tested on 24 Petri dishes. Most of the differences in the size of the clear zone exist because the samples themselves differ (some are stronger: A and C are the strongest, F the weakest). A smaller part exists because dishes differ slightly from one another, and very little is left to random measurement error. The test is therefore reliable: repeating a sample on a new dish gives a result close to the earlier one. But because only six samples were tried, we must be careful when generalising to all possible samples, and we must not treat the 144 numbers as 144 independent facts: they carry the information of about seven.

4.2 Limitations

General (d). What are the limitations of the model? (1) Only six samples: the sample variance is imprecise (CI 1.10 to 3.56) and the normality of sample effects cannot be verified; the estimate also depends strongly on sample F. (2) The model assumes additive effects; with one observation per cell, a plate by sample interaction is confounded with the residual. (3) One mildly extreme observation (retained after the sensitivity analysis). (4) Diameters are rounded to whole millimetres, so residuals are slightly discrete. (5) No plate-level covariates (agar batch, incubation position, date) were available. (6) Generalisation is limited to the conditions of this experiment.

5 Conclusion

A linear mixed model with crossed random intercepts for plate and sample fits the penicillin assay well (REML, no singularity, assumptions acceptable). The overall mean inhibition diameter is 22.97 mm (95% CI 21.27 to 24.68). Sample-to-sample differences dominate (78.5% of the variance), plate-to-plate variation is modest (15.1%) and residual error is small (6.4%). Ignoring the crossed structure understates the uncertainty of the mean by a factor of about five. If the six samples are the objects of interest, treating them as fixed gives precise pairwise comparisons (all pairs differ except A and C, and D and E); if they represent a population, the random formulation is the appropriate one. A larger number of samples and replicated cells would allow a plate by sample interaction to be estimated.

References

• Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1-48.
• Davies, O. L., and Goldsmith, P. L. (1972). Statistical Methods in Research and Production (4th ed.). Oliver and Boyd.
• Kuznetsova, A., Brockhoff, P. B., and Christensen, R. H. B. (2017). lmerTest package: Tests in linear mixed effects models. Journal of Statistical Software, 82(13), 1-26.
• Lüdecke, D., Ben-Shachar, M. S., Patil, I., Waggoner, P., and Makowski, D. (2021). performance: An R package for assessment, comparison and testing of statistical models. Journal of Open Source Software, 6(60), 3139.
• Nakagawa, S., and Schielzeth, H. (2013). A general and simple method for obtaining R² from generalized linear mixed-effects models. Methods in Ecology and Evolution, 4(2), 133-142.
• Nieuwenhuis, R., te Grotenhuis, M., and Pelzer, B. (2012). influence.ME: Tools for detecting influential data in mixed effects models. The R Journal, 4(2), 38-47.
• Pinheiro, J. C., and Bates, D. M. (2000). Mixed-Effects Models in S and S-PLUS. Springer.
• Self, S. G., and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605-610.
• Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A., and Smith, G. M. (2009). Mixed Effects Models and Extensions in Ecology with R. Springer.