What changes when the same sample is measured more than once?
Throughout this series, we have treated the 240 observations as observations of relative bias.
But there is an important feature of the experimental design that we have not yet modelled explicitly.
There are not 240 independent samples.
There are 60 environmental samples, and each sample was analysed with all four methods.
So the data have the structure:
\[ 60\ \text{samples} \times 4\ \text{methods} = 240\ \text{observations} \]
The four measurements belonging to the same sample share a common origin.
They may therefore be more similar to one another than measurements from different samples.
This changes the question.
What does it mean to compare methods when every method is measured on the same samples?
This is the domain of repeated measurements and blocked designs.
It is also the point where the assumption of independent observations becomes impossible to ignore.
Before fitting another model, we need to distinguish three concepts:
Our dataset contains 60 samples.
Each sample contributes four method-specific measurements.
The design can therefore be represented as:
\[ \begin{array}{c|cccc} \text{Sample} & \text{Reference} & \text{Method A} & \text{Method B} & \text{Method C}\\ \hline 1 & Y_{11} & Y_{12} & Y_{13} & Y_{14}\\ 2 & Y_{21} & Y_{22} & Y_{23} & Y_{24}\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 60 & Y_{60,1} & Y_{60,2} & Y_{60,3} & Y_{60,4} \end{array} \]
The important point is that the four values in each row belong to the same sample.
They are therefore paired within sample.
This is different from a design in which 240 unrelated samples were randomly divided among four methods.
Many of the models we have used so far can be written schematically as:
\[ Y_i = \mu + \text{systematic effects} + \varepsilon_i \]
with:
\[ \operatorname{Cov}(\varepsilon_i,\varepsilon_j)=0 \qquad i\neq j \]
The residuals are assumed to be independent.
But if two measurements come from the same sample, this may not be realistic.
Suppose Sample 17 has unusually high organic matter, an unusual mineral composition, or some other characteristic that makes its relative bias systematically different from the other samples.
That characteristic affects all four method measurements.
The residuals are therefore not simply four unrelated errors.
They share information about the same underlying sample.
The issue is not that we have “too many observations”.
We have 240 measurements.
The issue is that 240 measurements do not contain the same amount of independent information as 240 measurements from unrelated samples.
The first step is to make the repeated structure visible.
library(here)
library(data.table)
library(ggplot2)
dat <- fread(
here("data", "dataset_v01_clean_balanced.csv")
)
dat[, Method := factor(
Method,
levels = c("Reference", "Method_A", "Method_B", "Method_C")
)]
dat[, Matrix := factor(
Matrix,
levels = c("Fill", "Soil", "Sediment")
)]
dat[, .N, by = SampleID][, .N, by = N] N N
<int> <int>
1: 4 60There are 60 unique samples.
Every sample contributes exactly four observations.
This is the defining structure of the dataset for the next two chapters.
A useful way to see the structure is to connect measurements belonging to the same sample.
We do not need to display all 60 samples.
A random subset is easier to read.
set.seed(42)
sample_plot <- dat[
SampleID %in% sample(unique(SampleID), 12)
]
ggplot(
sample_plot,
aes(x = Method, y = RelativeBias_pct, group = SampleID)
) +
geom_line(alpha = 0.35) +
geom_point(size = 2) +
labs(
x = NULL,
y = "Relative Bias (%)",
title = "Repeated measurements within the same sample",
subtitle = "Each line connects the four methods measured on one sample"
) +
theme_minimal(base_size = 12)
The lines are not random noise: some samples tend to have relatively high bias across several methods, others tend to have relatively low bias.
The method effect is superimposed on this sample-to-sample variation.
This is exactly the structure that an independent-observations model fails to represent explicitly.
There is a simpler statistical idea that already captures part of this structure.
Suppose we wanted to compare Method B with Method C only.
Because both methods were applied to the same 60 samples, we could calculate:
\[ D_k = Y_{k,B} - Y_{k,C} \]
for each sample \(k\).
The comparison is then based on 60 within-sample differences rather than treating the two methods as unrelated groups.
This is the logic of a paired t-test.
wide_bc <- dcast(
dat[Method %in% c("Method_B", "Method_C")],
SampleID ~ Method,
value.var = "RelativeBias_pct"
)
t.test(
wide_bc$Method_B,
wide_bc$Method_C,
paired = TRUE
)
Paired t-test
data: wide_bc$Method_B and wide_bc$Method_C
t = 8.7003, df = 59, p-value = 3.62e-12
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
3.506292 5.600875
sample estimates:
mean difference
4.553583 The paired analysis asks a very specific question:
Within the same sample, is the expected bias different between Method B and Method C?
This is often more informative than an independent-samples comparison because sample-specific variation is removed from the comparison.
The general principle is important:
When measurements are naturally paired, the analysis should respect the pairing.
But a paired t-test is not enough for our full factorial design.
We have four methods, three matrices, and an interaction.
We need a model that can represent all of these structures simultaneously.
There is another way to describe the same design.
Each sample can be regarded as a block.
Within a block, all four methods are compared on the same material.
This has a major experimental advantage.
Suppose Sample 17 has unusually high organic matter.
Because all four methods are applied to Sample 17, that unusual characteristic affects all methods within the block.
The comparison between methods can therefore be made while holding the sample itself approximately constant.
This is the logic of a randomized block design: the sample acts as a blocking factor, the method is the treatment we want to compare.
One possible model is to include SampleID as a fixed effect:
\[ Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \gamma_k + \varepsilon_{ijk} \]
where:
In R:
fit_block <- lm(
RelativeBias_pct ~ Method * Matrix + SampleID,
data = dat
)
anova(fit_block)Analysis of Variance Table
Response: RelativeBias_pct
Df Sum Sq Mean Sq F value Pr(>F)
Method 3 13127.4 4375.8 528.8626 < 2.2e-16 ***
Matrix 2 804.8 402.4 48.6327 < 2.2e-16 ***
SampleID 57 818.6 14.4 1.7357 0.003603 **
Method:Matrix 6 185.9 31.0 3.7440 0.001587 **
Residuals 171 1414.9 8.3
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This model removes systematic differences between samples from the residual error.
It is therefore much closer to the actual experimental comparison than a model that treats SampleID as irrelevant.
But it also introduces 59 additional parameters.
That is not necessarily wrong.
It is simply not the most natural way to describe a population of samples.
We usually do not want to estimate a separate scientific parameter for Sample 1, Sample 2, Sample 3, and so on.
We want to describe how much samples vary from one another.
That is the motivation for random effects.
The distinction is conceptual.
A fixed effect is usually used when the levels themselves are the scientific objects of interest.
For example:
We want to compare these specific levels.
A random effect is used when the observed levels represent a sample from a broader population of possible levels, and our interest is primarily in the variability they introduce.
The 60 SampleIDs are a natural example.
We are not interested in whether Sample 17 is different from Sample 23.
We are interested in the fact that samples differ.
That is a variance question.
The same reasoning can apply to laboratory process variables.
Our dataset contains variables such as:
BatchDayAnalystInstrumentBut their presence in the dataset does not automatically mean that they should all become random effects.
The first question is:
Do observations sharing the same level of this variable plausibly share a source of dependence or common variation?
For example, measurements performed in the same analytical batch may share:
Measurements performed on the same day may also share common conditions.
These are plausible sources of clustering.
But there is another question:
Does the experimental design provide enough independent levels to estimate that variance component?
A random effect with only two or three levels is difficult to estimate reliably.
And if Day, Batch, Analyst and Instrument are strongly nested or confounded, adding all of them to the model may create an apparently sophisticated model that the design cannot actually support.
Therefore we will not add them simply because they are available.
They are candidates for random effects only if their experimental structure justifies them.
We will return to this issue in the next chapter.
The paired t-test is useful when the scientific question is a single comparison between two methods.
Our question is broader.
We want to estimate:
This is more naturally expressed as a hierarchical model.
The structure is:
\[ \text{measurement} \rightarrow \text{sample} \rightarrow \text{population of samples} \]
The mixed model provides a way to represent that hierarchy directly.
The word pseudoreplication is often used whenever repeated observations are analysed as if they were independent.
But the more precise issue in our dataset is first of all dependence.
The four measurements from a sample are not automatically invalid observations.
They are useful repeated measurements.
The problem arises when the analysis treats them as if they came from independent experimental units.
Whether this constitutes pseudoreplication depends on the experimental question and on what was independently randomized.
The practical lesson is simpler:
Before choosing a statistical test, identify the unit that was independently replicated or randomized.
A diagnostic plot cannot recover this information.
It comes from the experimental design.
Earlier we asked:
“Do methods differ?”
Then we asked:
“Does the method effect depend on matrix?”
Now the question becomes:
“Do methods differ within samples, while accounting for the fact that samples themselves vary?”
This is a more precise statement of the experimental comparison.
The statistical model should reflect it.
Repeated measurements are not automatically independent observations.
Our dataset contains 240 measurements, but they are organised into 60 samples with four measurements per sample.
That structure creates dependence.
A paired analysis can account for this when comparing two methods.
A blocked fixed-effects model can account for it by treating SampleID as a fixed blocking factor.
But when the scientific interest is in a population of samples rather than in the 60 observed samples individually, a random-effects formulation is more natural.
This leads directly to mixed-effects models.
The important lesson is not:
“Use a mixed model whenever you have repeated measurements.”
It is:
Identify the dependence structure created by the experimental design, then choose a model that represents the structure relevant to the scientific question.
We now know why SampleID matters.
The next question is how to represent sample-to-sample variation without estimating a separate fixed effect for every sample.
We will introduce a random intercept for SampleID.
Then we will ask:
“How much of the total variability comes from differences between samples, and how does acknowledging that structure change our inference about Method and Matrix?”