20.3 – Baseline correction

Introduction.

Signal processing is the analysis and manipulation of signals to extract meaningful information and improve data quality. signal processing is a crucial pre-processing step before analysis, as it involves cleaning and preparing raw data to improve its quality and highlight important features for accurate analysis. Common pre-processing tasks include filtering noise, filling in missing data, and feature extraction, all of which are done before the main analytical steps like feature extraction and classification can be performed.

A baseline refers to an initial measurement before an intervention. A starting point. Provides an objective comparison to judge whether or not an intervention has led to change.

Baseline drift is a gradual, slow shift in the signal’s zero-point over time, while intensity drift is a more general term for a change in a peak’s amplitude, often due to baseline changes. Noise is rapid, random fluctuation that obscures the signal itself. The key differences lie in their speed and effect: drift is a slow, low-frequency, long-term issue that shifts the entire baseline, whereas noise is a fast, high-frequency, short-term issue that adds random variation to the signal

For measures conducted over time, a baseline correction may be applied during initial data processing to correct for signal distortion, background noise, or baseline drift — the gradual change over time in what is expected to be measurement of an unchanging signal.

Many applications, numerous examples.

Quantitative PCR, qPCR: baseline correction is the process of identifying and subtracting the background fluorescence noise from the early cycles of a real-time PCR run to accurately measure the signal from specific DNA amplification (see Ruijter et al 2009).

Chromatography: baseline correction is a technique to remove background noise and drift from a chromatogram to make peaks more visible and accurate for analysis (see Niezen et al 2022).

Standard and basal metabolic rate by indirect calorimetry: baseline drift is an error where the instrument’s signal gradually changes over time, independent of the subject’s actual oxygen consumption (), while baseline correction is a data processing technique used to computationally adjust the recorded data to counteract this drift and improve accuracy (see Hayes et al 1992; Lighton 2017).

Colorimetric spectrometry: baseline drift is an unwanted phenomenon where the signal’s baseline gradually shifts over time or wavelength due to instrumental or environmental factors, while baseline correction is a data processing technique used to computationally remove this drift and other background noise from the measured spectrum.

Statistical considerations.

When processing a biological signal such as an electromyogram (EMG), it’s important to remember that even the baseline—the part of the recording where we assume “nothing is happening” — is only an estimate of the true baseline noise level. To account for this, choose two kinds of time windows: a baseline window (before the event of interest) and a signal window (during the activity you want to measure).

The goal is to compare these windows in a way that fairly adjusts for natural fluctuations in the baseline. One common approach is a regression-weighted correction, which simply means using a statistical line that represents how the baseline trends upward or downward over time, then adjusting the signal based on that line rather than assuming the baseline is perfectly flat. Another approach is to use a spline, which is a smooth, flexible curve that adapts to gradual changes in the baseline. Splines can correct for slow drifts in the recording without over-correcting the actual signal. Together, these methods help ensure that any “activity” you detect is more likely to be real muscle activation and not just shifts in the baseline.

Reporting detrended data ensures that your analysis reflects true physiological changes rather than technical artifacts like sensor drift or baseline wander. It establishes an unvarnished baseline that guarantees reproducibility across trials, allowing you to reliably compare a subject’s metabolic capacity over multiple tests regardless of shifting day-to-day instrument errors.

Example, EMG recording.

To illustrate, consider a myogram signal trace (EMG) recorded over several minutes. The trace will show fluctuations in electrical activity over time, which reflects muscle rest, contraction, and relaxation. The trace will exhibit a baseline at rest, spikes or bursts of activity during contractions, and varying levels of intensity depending on the muscle’s effort. An increase in the frequency and amplitude of these spikes indicates a stronger, more forceful contraction, while a period of no activity will show as a flat line. In a real-world scenario with an active muscle or a prolonged recording, the trace is typically corrupted by both noise (high-frequency, random variations) and baseline drift (slow, low-frequency shifts from the zero point). Among several statistics, analyst may calculate from the signal (1) the Root Mean Square (RMS) Amplitude, a measure of the signal’s magnitude and represents the overall intensity of muscle activity, (2) Mean Power Frequency (MPF), or the the frequency domain analysis of myogram signals. The MPF is calculated to assess muscle fatigue (fatigue often causes a shift to lower frequencies), and others (eg, Potvin and Bent 1997, Smilios et al 2010). Principle to the analysis, moving average or LOESS approaches may be used to smooth the data, reduce noise, and identify underlying trends or patterns in muscle activity over the multiple-minute duration.

To demonstrate pre-processing steps to correct for baseline drift, we need a data set. Here, we simulate a couple of myogram-like traces.

R code to simulate myogram data with baseline drift and random walk noise. A random walk tends to wander and works well for simulating biological drift.

# Example
# Define simulation parameters
duration_sec <- 5 # Total duration in seconds
sampling_rate_hz <- 1000 # Sampling rate (1000 Hz = 1 ms interval)
peak_center_time <- 2.5 # Time of the muscle contraction peak
amplitude <- 5 # Peak amplitude of the myogram
drift_magnitude <- 0.01 # Factor to control the magnitude of baseline drift
noise_level_sd <- 0.1 # Standard deviation of the white noise

Generate the data

my_data <- simulate_myogram_with_drift(
duration = duration_sec,
hz = sampling_rate_hz,
peak_time = peak_center_time,
peak_amplitude = amplitude,
drift_factor = drift_magnitude,
noise_sd = noise_level_sd
)
head(my_data)

Plot the simulated data

plot(my_dataSignal, type = 'l', col = 'blue',
main = "Simulated Myogram Data with Baseline Drift",
xlab = "Time", ylab = "Signal Value")
lines(my_dataDrift, col = 'red', lty = 2) # Overlay the drift line
legend("topleft", legend = c("Simulated Myogram", "Baseline Drift"),
col = c("blue", "red"), lty = c(1, 2), cex = 0.8)
# ylim = range(myogram_data, baseline_drift)

which gives us graph like Figure 1.

Figure 1. Simulated myogram data with baseline drift.

Alternatively, use ggplot (Fig 2).

library(ggplot2)
ggplot(my_data, aes(x = Time, y = Signal)) +
geom_line(color = "blue") +
geom_line(aes(y = Drift), color = "red", linetype = "dashed", alpha = 0.6) +
geom_line(aes(y = TrueSignal), color = "green", linetype = "dotted", alpha = 0.8) +
labs(title = "Simulated Myogram Data with Baseline Drift",
x = "Time (seconds)",
y = "Signal Amplitude") +
theme_minimal() +
scale_color_manual(values = c("blue", "red", "green"),
labels = c("Total Signal", "Baseline Drift", "True Myogram")) +
theme(legend.position = "bottom")

# You can access the data for further analysis using the 'my_data' data frame
head(my_data)
}

End R code.

version 2
# R code to simulate myogram data with baseline drift

# 1. Define simulation parameters
set.seed(123) # for reproducibility
n_samples <- 500 # number of data points (time steps)
time <- 1:n_samples
baseline_start <- 5
drift_rate <- 0.01 # slope of the linear drift
signal_amplitude <- 2
noise_sd <- 0.5 # standard deviation of the noise

# 2. Simulate the baseline drift
# A simple linear drift is used here. You could also use a random walk (RW).
baseline_drift <- baseline_start + drift_rate * time

# 3. Simulate the biological signal (myogram activity)
# Using a sinusoidal function as an example of a rhythmic signal
biological_signal <- signal_amplitude * sin(time * 0.1)

# 4. Simulate random noise (white noise)
noise <- rnorm(n_samples, mean = 0, sd = noise_sd)

# 5. Combine all components to get the final simulated myogram data
myogram_data <- baseline_drift + biological_signal + noise

# 6. Create a data frame for plotting and analysis
sim_data <- data.frame(Time = time, Signal = myogram_data)

# 7. Visualize the data using base R graphics or ggplot2
plot(sim_dataSignal, type = 'l', col = 'blue',
main = "Simulated Myogram Data with Baseline Drift",
xlab = "Time", ylab = "Signal Value", ylim = range(myogram_data, baseline_drift))
lines(sim_data\dot{V}O_{2}maxBMR$, values within those 20-second windows, you need to calculate the average of the raw signal values captured inside the logical masks.

Below is the updated base R code that detrends the data first, identifies the true physiological extremes, computes the window means, and prints them directly to the console alongside the plot.

R code.
# 1. Setup and Simulate Data with Linear Drift
set.seed(42)
n_points <- 1000
time_sec <- seq(0, 100, length.out = n_points)

signal <- sin(seq(0, 10, length.out = n_points)) + rnorm(n_points, sd = 0.1)
drift <- 0.01 * time_sec 
data_stream <- signal + drift

# 2. Detrend the data using linear regression residuals
drift_model <- lm(data_stream ~ time_sec)
detrended_data <- residuals(drift_model)

# 3. Identify indices based on the CLEAN, detrended data
zenith_idx <- which.max(detrended_data)
nadir_idx <- which.min(detrended_data)

zenith_time <- time_sec[zenith_idx]
nadir_time <- time_sec[nadir_idx]

# 4. Define 20-second Windows (+/- 10s around index)
window_half_width <- 10

zenith_window_mask <- (time_sec >= (zenith_time - window_half_width)) & 
(time_sec <= (zenith_time + window_half_width))

nadir_window_mask <- (time_sec >= (nadir_time - window_half_width)) & 
(time_sec <= (nadir_time + window_half_width))

# Compute the mean of the actual signal values within these 20-second blocks
mean_zenith_val <- mean(data_stream[zenith_window_mask])
mean_nadir_val <- mean(data_stream[nadir_window_mask])

# Print metrics to the R console
cat("--- Physiological Window Metrics ---\n")
cat(sprintf("Zenith Window Center Time: %.2f seconds\n", zenith_time))
cat(sprintf("Mean Zenith Window Value: %.4f\n\n", mean_zenith_val))
cat(sprintf("Nadir Window Center Time: %.2f seconds\n", nadir_time))
cat(sprintf("Mean Nadir Window Value: %.4f\n", mean_nadir_val))
cat("------------------------------------\n")

# Expand right margin (5th value) to make clean room for the legend outside the axis box
old_pars <- par(no.readonly = TRUE)
par(mar = c(5, 4, 4, 8) + 0.1, xpd = TRUE)

plot(time_sec, detrended_data, type = "l", col = "purple", lwd = 1.5,
xlab = "Time (seconds)", ylab = "Detrended Residuals",
main = "Clean Detrended Signal & Target Window Averages")

# Shading the windows on the detrended baseline
rect(zenith_time - window_half_width, min(detrended_data),
zenith_time + window_half_width, max(detrended_data),
col = rgb(1, 0, 0, 0.12), border = NA)

rect(nadir_time - window_half_width, min(detrended_data),
nadir_time + window_half_width, max(detrended_data),
col = rgb(0, 0, 1, 0.12), border = NA)

# Highlight window segments
lines(time_sec[zenith_window_mask], detrended_data[zenith_window_mask], col = "darkred", lwd = 2.5)
lines(time_sec[nadir_window_mask], detrended_data[nadir_window_mask], col = "darkblue", lwd = 2.5)

# Plot exact peaks identified after detrending
points(zenith_time, detrended_data[zenith_idx], col = "red", pch = 19, cex = 1.8)
points(nadir_time, detrended_data[nadir_idx], col = "blue", pch = 19, cex = 1.8)

# Overlay horizontal dashed lines tracking detrended window means
abline(h = mean_zenith_detrended, col = "darkred", lty = 2, lwd = 1.5)
abline(h = mean_nadir_detrended, col = "darkblue", lty = 2, lwd = 1.5)

# Position a clean, compact legend cleanly outside the plot area on the right margin
legend("topleft", inset = c(1.02, 0.2), bty = "n",
lty = c(1, 1, 1, 1, 2, 2), pch = c(NA, NA, NA, 19, NA, NA),
lwd = c(1.5, 2.5, 2.5, NA, 1.5, 1.5), 
col = c("purple", "darkred", "darkblue", "black", "darkred", "darkblue"),
legend = c("Detrended Signal", 
"Zenith Window", 
"Nadir Window", 
"Peak Points",
sprintf("Mean Zen: %.2f", mean_zenith_detrended),
sprintf("Mean Nad: %.2f", mean_nadir_detrended)))

# Restore original plotting space settings
par(old_pars)
End R code.
R output.
The mean Zenith Window Value (Raw) was 1.6102; after accounting for the baseline drift, the mean detrended Zenith was 0.6605 -- without correction we would report a value over 200% higher than the degtrended value. Similarly, the mean Nadir Window Value (Raw) was -0.3690; after detrending, mean degtrended Nadir was -1.0276, a value 36% lower than the uncorrected value.

Figure 4. Detrended zenith (red) and nadir (blue) highlighted for our simulated physiolical trace.
Questions

Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help

References
Hayes, J. P., Speakman, J. R., & Racey, P. A. (1992). Sampling bias in respirometry. Physiological Zoology, 65, 604–619.
Lighton, J. R. B. (2017). Limitations and requirements for measuring metabolic rates: A mini review. European Journal of Clinical Nutrition, 71(3), 301–305.
Liland, K. H. (2015). 4S Peak Filling–baseline estimation by iterative mean suppression. MethodsX, 2, 135-140.
Liland, K. H., & Mevik, T. A. B. H. (2011). Optimal baseline correction for multivariate calibration using open-source software. Life Science Instruments, (3), 7.
Niezen, L. E., Schoenmakers, P. J., & Pirok, B. W. J. (2022). Critical comparison of background correction algorithms used in chromatography. Analytica Chimica Acta, 1201, 339605.
Ruijter, J. M., Ramakers, C., Hoogaars, W. M. H., Karlen, Y., Bakker, O., van den Hoff, M. J. B., & Moorman, A. F. M. (2009). Amplification efficiency: Linking baseline and bias in the analysis of quantitative PCR data. Nucleic Acids Research, 37(6), e45.
Potvin, J. R., & Bent, L. R. (1997). A validation of techniques using surface EMG signals from dynamic contractions to quantify muscle fatigue during repetitive tasks. Journal of Electromyography and Kinesiology, 7(2), 131–139.
Smilios, I., Hakkinen, K., & Tokmakidis, S. P. (2010). Power Output and Electromyographic Activity During and After a Moderate Load Muscular Endurance Session. The Journal of Strength and Conditioning Research, 24(8), 2122–2131.

Chapter 20 contents

Additional topics
Area under the curve
Peak detection
Baseline correction
Surveys
Time series
Dimensional analysis
Estimating population size
Diversity indexes
Survival analysis
Growth equations and dose response calculations
Plot a Newick tree
Phylogenetically independent contrasts
How to get the distances from a distance tree
Binary classification
Meta-analysis

/MD

				

12.3 – Fixed effects, random effects, and agreement

Introduction

Within discussions of one-way ANOVA models the distinction between two general classes of models needs to be made clear by the researcher. The distinction lies in how the levels of the factor are selected. If the researcher selects the levels, then the model is a Fixed Effects Model, also called a Model I ANOVA. On the other hand, if the levels of the factor were selected by random sampling from all possible levels of the factor, then the model is a Random Effects Model, also called a Model II ANOVA. This page is about such models and I’ll introduce the intraclass correlation coefficient, abbreviated as ICC, as a way to illustrate applications of Model II ANOVA.

Here’s an example to help the distinction. Consider an experiment to see if over the counter painkillers are as good as prescription pain relievers at reducing numbers of migraines over a six week period. The researcher selects Tylenol^®, Advil^®, Bayer^® Aspirin, and Sumatriptan (Imitrex^®), the latter an example of a medicine only available by prescription. This is clearly an example of fixed effects; the researcher selected the particular medicines for use.

Random effects, in contrast, implies that the researcher draws up a list of all over the counter pain relievers and draws at random three medicines; the researcher would also randomly select from a list of all available prescription medicines.

Fixed effects are probably the more common experimental approaches. To be complete, there is a third class of ANOVA called a Mixed Model or Model III ANOVA, but this type of model only applies to multidimensional ANOVA (eg, two-way ANOVA or higher), and we reserve our discussion of the Model III until we discuss multidimensional ANOVA (Table 1).

Table 1. ANOVA models

Although the calculations for the one-way ANOVA under Model I or Model II are the same, the interpretation of the statistical significance is different between the two.

In Model I ANOVA, any statistical difference applies to the differences among the levels selected, but cannot be generalized back to the population. In contrast, statistical significance of the Factor variable in Model II ANOVA cannot be interpreted as specific differences among the levels of the treatment factor, but instead, apply to the population of levels of the factor. In short, Model I ANOVA results apply only to the study, whereas Model II ANOVA results may be interpreted as general effects, applicable to the population.

This distinction between fixed effects and random effects can be confusing, but it has broad implications for how we interpret our results in the short-term. This conceptual distinction between how the levels of the factor are selected also has general implications for our ability to acquire generalizable knowledge by meta-analysis techniques (Hunter and Schmidt 2000). Often we wish to generalize our results: we can do so only if the levels of the factor were randomly selected in the first place from all possible levels of the factor. In reality, this may not often be the case. It is not difficult to find examples in the published literature in which the experimental design is clearly fixed effects (ie, the researcher selected the treatment levels for a reason), and yet in the discussion of the statistical results, the researcher will lapse into generalizations.

Random Effects Models and agreement

Model II ANOVA is common in settings in which individuals are measured more than once. For example, in behavioral science or in sports science, subjects are typically measured for the response variable more than once over a course of several trials. Another common setting of Model II ANOVA is where more than one raters are judging an event or even a science project. In all of these cases what we are asking is about whether or not the subjects are consistent, in other words, we are asking about the precision of the instrument or measure.

In the assessment of learning by students, for example, different approaches may be tried and the instructor may wish to investigate whether the interventions can explain changes in test scores. To assess instrument validity — an item score from a grading rubric (Goertzen and Klaus 2023), tumor budding as biomarker for colorectal cancer (Bokhorst et al 2020) — agreement or concordance among two or more raters or judges (inter-rater reliability) needs to be established (same citations and references therein). For binomial data, we would be tempted to use the phi coefficient, Is it actually reliable Cohen’s kappa (two judges) or Fleiss Kappa (two or more judges) is a better tool. Phi coefficient and both assessment scores were introduced in Chapter 9.2.

For ratio scale data, agreement among scores of  two or more scores from judges There is an enormous number of articles on reliability measures in the social sciences and you should be aware of a classical paper on reliability by Shrout and Fleiss (1979) (see also McGraw and Wong, 1996). Both the ICC and the Pearson product moment correlation, r, which we will introduce in Chapter 16, are measures of strength of linear association between two ratio scale variables (Jinyuan et al 2016). But ICC is more appropriate for association between repeat measures of the same thing, eg, repeat measures of running speed or agreement between judges. In contrast, the product moment correlation can be used to describe association between any two variables, eg, between repeat measures of running speed, but also between say running speed and maximum jumping height. The concept of repeatability of individual behavior or other characteristics is also a common theme in genetics and so you should not be surprised to learn that the concept actually traces to RA Fisher and his invention of ANOVA. Just as in the case of the sociology literature, there are many papers on the use and interpretation of repeatability in the evolutionary biology literature (eg, Lessels and Boag 1987; Boake 1989; Falconer and Mackay 1996; Dohm 2002; Wolak et al 2012).

Agreement and repeatability are Model II (random effects) because the subjects, items, or raters being studied are treated as a random sample from a larger population rather than specific, fixed levels of interest. In reliability studies, subjects are not chosen for their specific characteristics but as representatives of the population variability. The goal is to infer whether a measurement tool works for any future subject or rater, not just the ones tested in the study. These factors are designed to estimate variance components and generalize results beyond the specific sample used. Put another way, the goal of Model II is often to estimate how much variance is due to differences between subjects versus within subjects (repeatability), which requires treating the subject effect as random.

There are many ways to analyze these kinds of data but a good way is to treat this problem as a one-way ANOVA with Random Effects. Thus, the Random Effects model permits the partitioning of the variation in the study into two portions: the amount that is due to differences among the subjects or judges or intervention versus the amount that is due to variation within the subjects themselves. The Factor is the Subjects and the levels of the factor are how ever many subjects are measured twice or more for the response variable.

If the subjects performance is repeatable, then the Mean Square Between (Among) Subjects, , component will be greater than the Mean Square Error component, , of the model. There are many measures of repeatability or reliability, but the intraclass correlation coefficient or ICC is one of the most common. The ICC may be calculated from the Mean Squares gathered from a Random Effects one-way ANOVA. ICC can take any value between zero and one.

where

and

B and W refers to the among group (between or among groups mean square) and the within group components of variation (error mean square), respectively, from the ANOVA, MS refers to the Mean Squares, and k is the number of repeat measures for each experimental unit. In this formulation k is assumed to be the same for each subject.

By example, when a collection of sprinters run a race, if they ran it again, would the outcome be the same, or at least predictable? If the race is run over and over again and the runners cross the finish lines at different times each race, then much of the variation in performance times will be due to race differences largely independent of any performance abilities of the runners themselves and the Mean Square Error term will be large and the Between subjects Mean Square will be small. In contrast, if the race order is preserved race after race: Jenny is 1st, Ellen is 2nd, Michael is third, and so on, race after race, then differences in performance are largely due to individual differences. In this case, the Between-subjects Mean Square will be large as will the ICC, whereas the Mean Square for Error will be small, and so will the ICC.

Can the intraclass correlation be negative?

In theory, no. Values for ICC range between zero and one. The familiar Pearson product moment correlation, Chapter 16, takes any value between -1 and +1. However, in practice, negative values for ICC will result if .

In other words, if the within group variability is greater than the among group variability, then negative ICC are possible. The question is whether negative estimates are ever biologically relevant or, simply result of limitations of the estimation procedure. Small ICC values and few repeats increases the risk of negative ICC estimates. Thus, a negative ICC would be “simply a(n) “unfortunate” estimate (Liljequist et al 2019). However, there may be situations in which negative ICC is more than unfortunate.  For example, repeatability, a quantitative genetics concept (Falconer and Mackay 1996), can be negative (Nakagawa and Schielzeth 2010).

Repeatability reflects the consistency of individual traits over time (or spatially, eg, scales on on a snake), is defined as the ratio of within individual differences and among individuals differences. Thus, within individual variance is attributed to general environmental effects, , whereas variance among individuals is both genetic and environmental. Scenarios for which repeated measures on the same individual may occur under environmental conditions, termed specific environmental effects or by Falconer and Mackay (1996), that negatively affect performance are not hard to imagine — it remains to be determined how common these conditions are and whether or not estimates of repeatability are negative.

ICC Example

I extracted 15 data points from a figure about nitrogen metabolism in kidney patients following treatment with antibiotics (Figure 1, Mitch et al 1977). I used a web application called WebPlot Digitizer (https://apps.automeris.io/wpd/), but you can also accomplish this task within R via the digitize package. I was concerned about how steady my hand was using my laptop’s small touch screen, a problem that very much can be answered by thinking statistically, and taking advantage of the ICC. So, rather than taking just one estimate of each point I repeated the protocol for extracting the points from the figure three times generating a total of three points for each of the 15 data points or 45 points in all. How consistent was I?

Let’s look at the results just for three points, #1, 2, and 3.

In the R script window enter

points = c(1,2,3,1,2,3,1,2,3)

Change points to character so that the ANOVA command will treat the numbers as factor levels.

points = as.character(points)
extracts = c(2.0478, 12.2555, 16.0489, 2.0478, 11.9637, 16.0489, 2.0478, 12.2555, 16.0489)

Make a data frame, assign to an object, eg, “digitizer”

digitizer = data.frame(points, extracts)

The dataset “digitizer” should now be attached and available to you within Rcmdr. Select digitizer data set and proceed with the one-way ANOVA.

Output from oneway ANOVA command

Model.1 <- aov(extracts ~ points, data=digitizer)

 summary(Model.1)
            Df    Sum Sq   Mean Sq  F value     Pr(>F)    
points       2 313.38895 156.69448 16562.49 5.9395e-12 ***
Residuals    6   0.05676   0.00946                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

numSummary(extracts , groups = points, 
+   statistics=c("mean", "sd"))
         mean           sd data:n
1  2.04780000 0.0000000000      3
2 12.15823333 0.1684708085      3
3 16.04890000 0.0000000000      3

End R output. We need to calculate the ICC

I’d say that’s pretty repeatable and highly precise measurement!

But is it accurate? You should be able to disentangle accuracy from precision based on our previous discussion (Chapter 3.5), but now in the context of a practical way to quantify precision.

ICC calculations in R

We could continue to calculate the ICC by hand, but better to have a function. Here’s a crack at the function to calculate ICC along with a 95% confidence interval.

myICC <- function(m, k, dfN, dfD) {
 testMe <- anova(m) 
 MSB <- testMe$"Mean Sq"[1]
 MSE <- testMe$"Mean Sq"[2]
 varB <- MSB - MSE/k
 ICC <- varB/(varB+MSE)
 fval <- qf(c(.025), df1=dfN, df2=dfD, lower.tail=TRUE)
 CI = (k*MSE*fval)/(MSB+MSE*(k-1)*fval)
 LCIR = ICC-CI
 UCIR = ICC+CI
 myList = c(ICC, LCIR, UCIR)
 return(myList)
}

The user supplies the ANOVA model object (eg, Model.1 from our example), k, which is the number of repeats per unit, and degrees of freedom for the among groups comparison (2 in this example), and the error mean square (6 in this case). Our example, run the function

m2ICC = myICC(Model.1, 3, 2,6); m2ICC

and R returns

[1] 0.9999396 0.9999350 0.9999442

with the ICC reported first, 0.9999396, followed by the lower limit (0.9999350) and the upper limit (0.9999442) of the 95% confidence interval.

Instead of writing your own function, several R packages can calculate the intraclass correlation coefficient (ICC) and its variants. These packages include: irr, psy, psych, and rptR.  For complex experiments involving multiple predictor variables, these packages are helpful for obtaining the correct ICC calculation (cf Shrout and Fleiss 1979; McGraw and Wong 1996). For the one-way ANOVA it is easier to just extract the information you need from the ANOVA table and run the calculation directly. We do so for a couple of examples.

Example: Are marathon runners consistent more consistent than my commute times?

Figure 1.Honolulu Marathon 2024 participant. Image credit: Pdubs.94, licensed under CC BY 4.0 (Creative Commons Attribution 4.0 International License), via Wikipedia.

A marathon is 26 miles, 385 yards (42.195 kilometers) long (World Athletics.org). And yet, tens of thousands of people choose to run in these events. For many, running a marathon is a one-off, the culmination of a personal fitness goal. For others, it’s a passion and a few are simply extraordinary, elite runners who can complete the courses in 2 to 3 hours (Table 3). That’s about 12.5 miles per hour. For comparison, my 20 mile commute on the H1 freeway on Oahu typically takes about 40 minutes to complete, or 27 miles per hour (Table 2, yes, I keep track of my commute times, per Galton’s famous maxim: “Whenever you can, count”).

Figure 2. View west along Interstate H-1 (Lunalilo Freeway) at 6:30AM on a weekday, Honolulu, Oahu, Hawaii. Image used with permission, credit: S. Dohm.

Table 2. A sampling of average commute speeds, miles per hour (mph),  on the H1 freeway during DrD’s morning commute (distance about 19.3 miles).

Try on your own.

Calculate the ICC for my commute speeds.

Run the one-way ANOVA to get the necessary mean squares and input the values into our ICC function. We have

require(psych)
m2ICC = myICC(AnovaModel.1, 4, 4,11); m2ICC
[1] 0.7390535 0.6061784 0.8719286

Repeatability, as estimated by the ICC, was 0.74 (95% CI 0.606, 0.872), for repeat measures of commute times.

We can ask the same about marathon runners — how consistent from race to race are these runners? The following data are race times drawn from a sample of runners who completed the Honolulu Marathon in both 2016 and 2017 in 2 to 3 hours (times recorded in minutes). In other words, are elite runners consistent?

Table 3. Honolulu marathon running times (min) for eleven repeat, elite runners.

Stacked worksheet, P refers to person, therefore P1 is person ID 1, P2 is person ID 2, and so on.

After running a one-way ANOVA, here are results for the marathon runners

m2ICC = myICC(Model.1, 2, 10,11); m2ICC
[1] 0.9780464 0.9660059 0.9900868

Repeatability, as estimated by the ICC, was 0.98 (95% CI 0.966, 0.990), for repeat measures of marathon performance. Put more simply, knowing what a runner did in 2016 I would be able to predict their 2017 race performance with high confidence, 98%!

And now, we compare: the runners are more consistent!

Clearly this is an apples-to-oranges comparison, but it gives us a chance to think about how we might make such comparisons. The ICC will change because of differences among individuals. For example, if individuals are not variable, if everyone finishes the race at the same time, then the ICC will be low (or even zero).  Because there is no variance to be attributed to individual differences, the proportion of total variance explained by those differences will be being minimal.

Visualizing repeated measures.

Waterfall plot.

If just two measures, then a profile plot will do. We introduced profile plots in the paired t-test chapter (Ch10.3). For more than two repeated measures, a line graph will do (see question 9 below). Another graphic increasingly used is the waterfall plot (waterfall chart), particularly for before and after measures or other comparisons against a baseline. For example, it’s a natural question for our eleven marathon runners — how much did they improve on the second attempt of the same course? Figure 3 shows my effort to provide a visual.

Figure 3. Simple waterfall plot of race improvement for Table 3 data. Dashed horizontal line at zero.

R code

Data wrangling: Create a new variable, eg, race21, by subtracting the first time from the second time. Multiply by -1 so that improvement (faster time in second race) is positive. Finally, apply descending sort so that improvement is on the left and increasingly poorer times are on the right.

df <- data.frame(
ID = rep(paste0("P", 1:11), each = 2),
Time = c(179.9,192.0,
129.9,130.8,
128.5,129.6,
179.4,179.7,
174.3,181.7,
177.2,176.2,
169.0,173.4,
174.1,175.2,
175.1,174.2,
163.9,175.9,
179.3,179.8)
)

# Convert to wide format
df_wide <- reshape(df, idvar = "ID", timevar = "ID",
direction = "wide")

# Extract first and second times manually
time_matrix <- matrix(df$Time, ncol = 2, byrow = TRUE)
colnames(time_matrix) <- c("Time1", "Time2")

# Calculate improvement (second - first)
base21 <- time_matrix[,2] - time_matrix[,1]

# Sort for waterfall display
base21 <- sort(base21, decreasing = TRUE)

# Waterfall plot
barplot(base21, col="blue", border="blue", space=0.5,
ylim=c(-20,20),
ylab="Improvement from first race (min)")
abline(h=0, lty=2)

Interpretation.

From the graph (Fig 3) we see that few of the racers improved their times the second time around — a few had decreased times.

Spaghetti plots.

Spaghetti plots, which display individual trajectories over time or multiple measurements, is a specialized line plot effective for visualizing the consistency of performance. The vertical thickness or spread of a single subject’s line across measurements indicates their level of repeatability. A “flat” line indicates high repeatability, while a volatile, zig-zagging line indicates low repeatability. Spaghetti plots are useful to visualize growth, decline, or stability over time for a single performance measure.

Figure 4. A spaghetti plot of average commute speeds from Table 2 data. 

R code

# 1. Enter the data from Table 2.
commute <- data.frame(
Week = 1:4,
Monday = c(28.5, 25.8, 26.2, 23.3),
Tuesday = c(23.8, 22.4, 22.6, 26.9),
Wednesday = c(28.5, 29.3, 24.9, 31.3),
Thursday = c(30.2, 26.2, 24.2, 26.2),
Friday = c(26.9, 27.7, 34.3, 30.2)
)

head(commute)

# 2. Set the X-axis points
days <- 1:5
day_names <- c("Mon", "Tue", "Wed", "Thu", "Fri")

# 3. Create empty plot
# Find range of Y to ensure all data fits
ymin <- min(commute[, 2:6]) - 2
ymax <- max(commute[, 2:6]) + 2

plot(NA, xlim = c(1, 5), ylim = c(ymin, ymax),
xlab = "Day of Week", ylab = "Average Speed (mph)",
main = "Weekly Commute Speeds", xaxt = "n")

# Add X-axis labels
axis(1, at = 1:5, labels = day_names)

# 4. Add lines for each week (spaghetti)
colors <- c("red", "blue", "darkgreen", "purple")
for (i in 1:nrow(commute)) {
# Plot lines
lines(days, commute[i, 2:6], col = colors[i], lwd = 2, type = "b", pch = 19)
}

# 5. Add legend
legend("topleft", legend = paste("Week", 1:4), col = colors, lwd = 2, pch = 19)

Parallel coordinates plots.

Parallel coordinates plot allows each entity to be represented by a unique, continuous line traced across several vertical axes (each representing a different measurement or variable). High repeatability is indicated when a line remains roughly horizontal across all axes, while poor repeatability is visualized as lines that zig-zag or change direction between measurements.

Figure 5. A parallel coordinates plot of average commute speeds from Table 2 data.

Both spaghetti plots and parallel coordinate plots can be used to visualize consistency of performance over multiple measures, but they differ significantly in their handling of data, scale, and context. Spaghetti plots are best for tracking within-subject variability over time, but quickly become a mess — like pasta tossed and sticking to a wall — when subject count increases. Parallel coordinate plots are a multivariate generalization — each line represents an individual’s profile across different tests, eg, accuracy, speed, VO2max.

R code.

# 1. Create the data frame
# see above

# 2. Extract only the numerical commute data for plotting
commute_data <- commute[, 2:6]

# 3. Setup the plot area
# Define range for y-axis
y_range <- range(commute_data)
# Define number of days (5)
num_days <- ncol(commute_data)

# Create empty plot
plot(NA, xlim = c(1, num_days), ylim = y_range, 
xaxt = "n", xlab = "Day of Week", ylab = "Average Speed (mph)",
main = "Average Commute Speed: Weekly Comparison")

# Add x-axis labels
axis(1, at = 1:num_days, labels = colnames(commute_data))

# 4. Draw vertical lines for each day
for (i in 1:num_days) {
lines(c(i, i), y_range, col = "grey", lty = 2)
}

# 5. Add lines for each week (repeated measures)
colors <- c("red", "blue", "darkgreen", "orange")
for (i in 1:nrow(commute_data)) {
lines(1:num_days, commute_data[i, ], col = colors[i], lwd = 2, type = "b", pch = 19)
}

# 6. Add legend
legend("topleft", legend = paste("Week", commute$Week), col = colors, lwd = 2, pch = 19)

Heat map.

Heat map quickly identify patterns, such as an individual’s performance improving, degrading, or remaining stable over time. For example, uniform, solid colors across rows (individuals) or columns (time/trials) indicate high consistency, while varied colors indicate poor repeatability.

Figure 6. A heat map of the commute speeds data set.

R code.

# 1. Create the data
# Week is excluded from matrix to only include numerical commute times
commute_data <- matrix(c(
28.5, 23.8, 28.5, 30.2, 26.9,
25.8, 22.4, 29.3, 26.2, 27.7,
26.2, 22.6, 24.9, 24.2, 34.3,
23.3, 26.9, 31.3, 26.2, 30.2
), nrow = 4, byrow = TRUE)

# 2. Add row and column names
rownames(commute_data) <- c("Week 1", "Week 2", "Week 3", "Week 4")
colnames(commute_data) <- c("Mon", "Tue", "Wed", "Thu", "Fri")

# 3. Create the heatmap
# Use t() to transpose so weeks are rows and days are columns
# Use rev() on rows to have Week 1 at the top
image_data <- t(apply(commute_data, 2, rev))

# Set up the plotting area (margins for labels)
par(mar = c(5, 6, 4, 2))

# Plot the matrix
# Option 1
image(1:nrow(image_data), 1:ncol(image_data), image_data,
col = heat.colors(12),
axes = FALSE, xlab = "Days", ylab = "",
main = "Commute Times Heatmap")

# Add axes labels
axis(1, at = 1:5, labels = colnames(commute_data))
axis(2, at = 1:4, labels = rev(rownames(commute_data)), las = 1)
box()

Interpretation.

The graphs (Fig 4 – 6) show a lack of trends in consistency of drive speeds. The heat map suggests more consistency on Monday — low speeds (Fig 4 & Fig 5), but that’s a stretch.

An example for you to work, from our Measurement Day

If you recall, we had you calculate length and width measures on shells from samples of gastropod and bivalve species. In the table are repeated measures of shell length, by caliper in mm, for a sample of Conus shells (Fig 7 and Table 4).

Figure 7. Conus shells, image by M. Dohm.

Repeated measures of the same 12 Conus shells.

Table 4. Unstacked dataset of repeated length measures on 12 shells

Maximum length from from anterior end to apex with calipers. Repeated measures were conducted blind to shell id, sampling was random over three separate time periods.

Questions

1. Consider data in Table 2, Table 3, and Table 4. True or False: The arithmetic mean is an appropriate measure of central tendency. Explain your answer.
2. Enter the shell data into R; Best to copy and stack the data in your spreadsheet, then import into R or R Commander. Once imported, don’t forget to change Sample to character, otherwise R will treat Sample as ratio scale data type. Run your one-way ANOVA and calculate the intraclass correlation (ICC) for the dataset. Is the shell length measure repeatable?
3. The marathon results in Table 2 are paired data, results of two races run by the same individuals. Make a profile plot for the data (see Ch10.3 for help).
4. A profile plot is used to show paired data; for two or more repeat measures, use a line graph instead.
   A. Produce a line graph for data in Table 4. The line graph is an option in R Commander (Rcmdr: Graphs Line graph…), and can be generated with the following simple R code:

lineplot(Sample, Measure1, Measure2, Measure3)

B. What trends, if any, do you see? Report your observations.
C. Make another line plot for data in Table 2. Report any trend and include your observations.

Quiz Chapter 12.3

Fixed effects, random effects, and agreement

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Chapter 12 contents

• Introduction
• The need for ANOVA
• One way ANOVA
• Fixed effects, random effects, and agreement
• ANOVA from “sufficient statistics”
• Effect size for ANOVA
• ANOVA posthoc tests
• Many tests one model
• References and suggested readings