Now hiring We're looking for a Social Media Coordinator — a paid, part-time position (~4 hours/month). Learn more & apply →

Lab Guide › Research Workflow

Data analysis

The statistical methods we use, what is required of every analysis, and where to learn each one.

Once your experiment finishes or your code converges, now comes the analysis!

This shouldn’t be the first time you’re thinking about this though as good experimental design involves a consideration of the analysis goals right from the outset.

This page is about how to choose a model, fit it, check it, and report what it supports (and what it doesn’t). How well you know any single number is covered in Error and uncertainty analysis.

The five stages

  1. Clean. Remove blanks, false triggers and artifacts. This process must be documented and procedural. Never edit raw data; save cleaned data separately (see the source-file rule in Data management).
  2. Explore. Make many ugly, unstyled x–y plots. One to two weeks. Work through the eight checks below.
  3. Model. The rest of this page.
  4. Quantify uncertainty. Error and uncertainty analysis.
  5. Plot. Figures.

Stage 2 is not stage 3. Exploring is how you find out what the data look like. Fitting is how you test what you claimed before you looked. Choosing your model after the exploratory plots is p-hacking, even when it does not feel like it.

What to check before you fit

The eight checks below are the data-exploration protocol from Zuur, Ieno & Elphick (2010). Run them in order, on every dataset, before you fit anything.

Formulate the hypothesis · run the experiment · collect the data

Data exploration

  1. Outliers in Y and Xboxplot · Cleveland dotplot
  2. Homogeneity of Yconditional boxplot
  3. Normality of Yhistogram · Q–Q plot
  4. Zero trouble in Yfrequency plot · corrgram
  5. Collinearity in XVIF · scatterplots · PCA
  6. Relationships between Y and Xmulti-panel scatterplots
  7. Interactionscoplots · conditional boxplots
  8. Independence of YACF · variogram · Y against time or space

Apply the statistical model

The data exploration protocol, after Zuur, Ieno & Elphick (2010), fig. 1. Step 8 is the one that decides your random-effect structure below.

Two of these do double duty on this page. Check 8, independence, is how you discover the non-independence that requirement 3 makes you model. Check 5, collinearity, is a common reason to reach for PCA rather than throwing every correlated predictor into one fit.

The requirements

Every analysis in this lab meets all seven. They are not negotiable and they are not ranked.

  1. The model is written down before the data are collected.
  2. The replicate unit is stated in words.
  3. Non-independence is in the model, not in a footnote.
  4. Residuals are checked, and the check is in the repository.
  5. The estimate and its interval are reported, not only the verdict.
  6. The report contains enough to refit the model.
  7. Uncertainty and inference are reported separately.

1 · Write the model down first

Before the first run, put a plain-text analysis plan in the project repository: response variable, predictors, random effects, and what result would count as support for each hypothesis. One paragraph is enough.

This is the cheapest protection that exists against fitting until something is significant (i.e., p-hacking). Makin & Orban de Xivry’s (2019) lists ten of the most common statistical mistakes, which is the shortest useful thing anyone in the lab can read on this.

2 · State the replicate unit

Write the sentence out: “n = 1 hawk, 5 trials per condition, ~4000 samples per trial; inference is about this individual.” Then check that the model’s degrees of freedom agree with it.

Counting correlated measurements as independent replicates is pseudoreplication (Hurlbert 1984; Lazic 2010). It remains the most common serious error in animal biomechanics, and a 2025 review of two decades of animal studies found it present in the majority of papers and increasing over time, despite statistical reporting improving across the same period. Better reporting does not fix this one. Naming the unit does.

3 · Put the non-independence in the model

Once the replicate unit is named, the structure follows. Repeated trials on one animal get a random intercept for trial. Repeated configurations of one wing specimen get a random intercept for the specimen. Species share ancestry, so cross-species comparisons carry the phylogeny.

Harrison et al. (2018) is the lab’s default reference for how to specify these models, and Bolker’s GLMM FAQ is where to go when one misbehaves.

4 · Check the residuals, and commit the check

Run the diagnostic, look at it, and leave the code in the repository so the next person can rerun it. Use DHARMa, which simulates scaled residuals for mixed models and tests dispersion and residual temporal, spatial and phylogenetic autocorrelation. A quantile–quantile plot you looked at once and did not save is not a check.

5 · Report the estimate, not only the verdict

A \(p\)-value alone says whether an effect is distinguishable from zero. It does not say how large it is, and that is usually the actual question. Report the estimate with a 95% interval alongside every test, and an effect size where one is meaningful.

This is the settled position across the field, not a stylistic preference: Wasserstein, Schirm & Lazar (2019) for the American Statistical Association, and Amrhein, Greenland & McShane (2019) in Nature, co-signed by over 800 researchers. Both argue for leading with magnitude and precision.

We keep \(\alpha = 0.05\) and two-sided tests as a convention for deciding what to discuss. We do not let it decide what to report.

6 · Report enough to refit

Model formula including the random terms, package and version, the degrees-of-freedom method, exact \(F\) and \(p\) (not “\(p < 0.05\)”), and the seed for anything resampled. See What to report.

7 · Keep uncertainty separate

Measurement uncertainty belongs in the error bars and the text, following Error and uncertainty analysis. Model uncertainty belongs in the confidence intervals. Do not let one stand in for the other. Where they genuinely have to combine, see Sharpening.

Choosing a model

Find the row that describes your data. Start there; justify anything more complicated.

Your data Start with Where it is explained
One value per condition, conditions independent lm() Schluter, Fit model → Simple linear regression
Repeated measurements on the same bird, specimen or model lmer(y ~ x + (1 \| ID)) Schluter, Fit model → Mixed models; Coding Club
Several individuals, generalising to the species lmer(), individual as random intercept Harrison et al. 2018
Multiple species gls(..., correlation = corPagel()), or PGLMM in MCMCglmm Schluter, Phylogenetic comparison; Hadfield’s course notes
Counts, proportions, binary outcomes glmer() or glmmTMB() Bolker’s GLMM FAQ
Two candidate functional forms (e.g. linear vs quadratic) Fit both, compare AIC — see warning below Harrison et al. 2018, §multi-model inference
Real uncertainty on both axes Major-axis or standardised major-axis regression (smatr) Schluter, Fit model → Correct for body size
A deterministic sweep with no sampling noise No test. Run a sensitivity study instead Harvey 2024
Many correlated shape or morphology variables, no single response Principal component analysis Schluter, Multivariate methods
A time series of coordinates with rhythmic structure Dynamic mode decomposition BirdDMD

Keep the candidate set small and justified. Every model in the set should correspond to a hypothesis you can state. Enumerating fifty polynomial combinations and taking the lowest AIC is a search, not an inference, and it will find structure in noise.

When the set is a genuine comparison, report Akaike weights as well as \(\Delta\)AIC: they say how much better, not just which.

What to report

Every fitted model in a paper, thesis chapter or committee meeting reports all of these. For how to word the result once you have the numbers, see Writing papers.

Item Example
Full model formula, random terms included avg_lift ~ state + perch_height + (1 \| trial_id)
Replicate structure in words “1 individual, 5 trials per condition”
Software and package versions R 4.4.1; lmerTest 3.1-3
Estimation and df method REML, Kenward–Roger
Test statistic with both degrees of freedom \(F_{1,17.06} = 42.43\)
Exact \(p\) \(p = 0.002\), or “\(p < 0.001\)” only below that
Estimate with a 95% interval \(-0.14\) [\(-0.21\), \(-0.07\)]
Effect size, where meaningful partial \(\eta^2 = 0.51\)
Model fit marginal and conditional \(R^2\)
Seed, for anything resampled set.seed(42)
Worked example: repeated measures on one animal

One bird, two feather conditions, three perch heights, several trials per combination. The design behind Martínez-Carmena et al. (2026), written out as a template you can adapt.

library(lmerTest)     # lmer with Kenward-Roger tests
library(DHARMa)       # residual diagnostics
library(performance)  # marginal / conditional R2
library(emmeans)      # marginal means

# 1. Fit. Fixed effects are the two designed factors;
#    the random intercept carries the repeated measures.
model_A <- lmer(avg_lift ~ factor_state + factor_perch + (1 | trial_id), data = data)

# 2. Check before you read anything off it.
res <- simulateResiduals(model_A)
plot(res)

# 3. Test. Kenward-Roger, because the sample is small.
anova(model_A, ddf = "Kenward-Roger")

# 4. Fit quality, both flavours.
r2(model_A)

# 5. Magnitudes and direction, on the response scale.
emmeans(model_A, ~ factor_state)

# 6. Intervals on the fixed effects.
confint(model_A, method = "profile")

Steps 4–6 are what turn a verdict into a result. They are the reason this analysis satisfies requirement 5 without any extra work.

Why Kenward–Roger. Luke (2017) showed that REML with \(F\)-tests using Kenward–Roger or Satterthwaite degrees of freedom gives the most accurate inference for fixed effects in linear mixed models. lmerTest defaults to Satterthwaite; Kenward–Roger is the more conservative choice at small \(n\) and is what we use.

Decomposition: finding structure in many variables

Sometimes there is no single response variable. A wing outline is fifty correlated coordinates; a motion-capture trial is twenty-four coordinates evolving in time. Decomposition methods turn that into a handful of axes or modes you can plot, interpret, and put into the models above.

The two we use, and how to tell them apart.

PCA finds the directions of greatest variance. It has no concept of time: shuffle the rows and you get the same answer.

DMD finds modes that each carry one frequency and one growth rate. Shuffle the rows and it is meaningless.

If your question contains “how much does shape vary, and along what axes”, you want PCA. If it contains “how fast”, “in what order”, or “what would the next wingbeat look like”, PCA cannot answer it and DMD can.

Principal component analysis

Schluter’s Multivariate methods page has the mechanics: prcomp(), scaling, scree plots, biplots, loadings and scores. Read that first. What follows is the part specific to how we use it.

1 · Get the variables onto a common scale before you start. prcomp() defaults to the covariance matrix (scale. = FALSE), which is right when everything shares units and has a comparable variance: usually after log-transforming a set of lengths. Mixing linear, area and mass measurements without correcting first will let whichever variable happens to have the largest numbers dominate PC1. Use scale. = TRUE (the correlation matrix) when the units genuinely differ. Schluter’s Preparing variables section gives the rule, including dividing logged areas by 2 and volumes by 3.

2 · If the variables are landmark coordinates, superimpose first. Raw digitised coordinates carry position, orientation and size, which will otherwise appear as the first few components. A generalised Procrustes superimposition removes all three and leaves shape. Use geomorph::gpagen(); it also handles sliding semi-landmarks along curves, which is what you want for a wing outline or a trailing edge.

3 · Decide how many components to keep, and decide it by a stated rule. Cumulative variance threshold, a break in the scree plot, whatever you like — but write the rule down before you look at the plot, and report it.

4 · Interpret by drawing, not by reading coefficients. Reconstruct and plot the shape at plus and minus a few standard deviations along each retained axis. A loadings table tells you almost nothing about what a wing is doing; two overlaid outlines tell you immediately.

5 · Using scores downstream is fine, with one caveat. PC scores make good response or predictor variables in the models above. But the model treats them as if they were measured, when they were estimated — the uncertainty in the components themselves is not carried through. Say so in the methods.

A small percentage is not automatically a small effect. In Gamble et al. (2020), PCA on 57 landmarks along a compliant trailing edge gave PC1 = 99.7% of shape change (driven by Reynolds number) and PC2 = 0.3% (driven by angle of attack). The second axis was 0.3% of the variance and still a real, interpretable, significant effect. Percentage of variance ranks the axes. It does not rank their importance to your question.

Report: how many variables and how many specimens went in, whether you used the covariance or correlation matrix, whether landmarks were Procrustes-aligned, variance explained by each retained component, the retention rule, and the sign convention.

Dynamic mode decomposition

DMD approximates a time series as a sum of modes, each with its own spatial shape, frequency and growth or decay rate:

\[\mathbf{x}(t) = \sum_{k=1}^{r} b_k\,\boldsymbol{\varphi}_k\, e^{\omega_k t}\]

That gives you something PCA cannot: a generative model. Because each mode carries a frequency, you can run it forward, synthesise a wingbeat that was never recorded, or reconstruct a trial from a chosen subset of modes and see what each one contributes.

Where to start. Our collaborator Lydia France has written BirdDMD, which wraps PyDMD with defaults chosen for biological motion-capture data. Work through the notebook gallery in order: it starts from a synthetic sum of two sinusoids where you can check DMD recovers the frequencies you put in, then moves to hawk flapping modes, turning manoeuvres, reconstruction error, and a generative model. The underlying study is France, Lapo & Kutz (2026), which decomposes flapping, turning, landing and gliding into a shared, low-order set of modes.

The knobs, and why they matter.

Setting What it does
n_modes (rank \(r\)) How many modes to keep. Too few and you lose real dynamics; too many and you fit noise
d (time-delay embedding) Stacks lagged copies of the data. Needed when you have fewer spatial channels than dynamics
eig_constraints conjugate_pairs forces eigenvalues into complex-conjugate pairs, which keeps the reconstruction real-valued and the modes interpretable as oscillations

Reconstruction error on the data you fitted proves nothing. Enough modes will reproduce any record. The number that means something is the error on a segment or a trial the fit never saw. Report that one.

Report: number of modes, delay embedding, any eigenvalue constraints, the frequency and growth rate of each retained mode, and held-out reconstruction error.

What DMD is and is not. It is a linear fit in a coordinate system chosen from the data. That is a strength when physics-based models rest on assumptions that do not survive real flight, and it is the limitation when you want a mechanism rather than a description. A DMD mode tells you the motion contains a coherent oscillation at that frequency and shape. It does not tell you why.

Sharpening the analysis

Everything above is required. Everything below is worth reaching for, and none of it is a prerequisite for a good paper. Pick what your question needs.

Say what you could have detected. With one or two animals, “power” is the wrong output. Simulate from a pilot fit (simr) and report the smallest effect the design could resolve. Do it before you collect, and the limitation becomes a stated scope rather than a caveat. Schluter’s Planning tools page covers the simpler pwr and power.t.test cases.

Answer the comparative question directly. emmeans and marginaleffects give contrasts on the response scale with intervals, and handle multiplicity adjustment when you are making many comparisons. Far more informative than a table of \(p\)-values.

Declare which tests are exploratory. When one study fits many models across many response variables at \(\alpha = 0.05\), some will clear the bar by chance. Either pre-specify a primary response, adjust, or label the rest exploratory in the text. Any of the three is defensible; silence is not.

Plot the fit, not just the data. visreg shows partial residuals against one predictor with the others accounted for, so a multi-predictor fit becomes something you can actually look at. It works with lmer fits.

Push measurement error through to the conclusion. This is the one place the two pages meet. Resample each point within its uncertainty range, refit, and check the conclusion still holds. Harvey et al. (2022) did this with 5,000 bootstrap draws over the centre-of-gravity error range, refitting the evolutionary model each time. Schluter’s Resample, bootstrap page has the mechanics, including boot and BCa intervals.

Show you had the power to choose the model. Selecting by AIC assumes the candidates are distinguishable at your sample size. Simulating under each and comparing likelihood-ratio distributions demonstrates it. Harvey et al. (2022) used pmc this way for the Ornstein–Uhlenbeck versus Brownian motion choice.

Consider a Bayesian fit when \(n\) is tiny. At two or three replicates an \(F\)-test is working with almost no degrees of freedom. A weakly informative prior gives a better-behaved interval. brms (Coding Club tutorial) or MCMCglmm. Optional, and never a way to rescue an analysis that failed frequentist assumptions.

Where to learn each method

Start with Dolph. The rest fill specific gaps he does not cover.

Dolph Schluter’s R tips pages — the lab’s default reference. Written for biologists, worked examples throughout.

If you need Go to
Linear models, ANOVA, mixed models, ML vs REML, singular fits Fit model → Read me, then Mixed models
Marginal means, post-hoc tests, visualising fits Fit model → Estimate magnitudes of effect
Body-size correction, errors in x, MA and SMA regression Fit model → Correct for body size
Simulating a design; power and sample size Planning tools
Bootstrap standard errors, BCa intervals, permutation tests Resample, bootstrap
Independent contrasts, phylogenetic GLS, Pagel’s \(\lambda\), OU Phylogenetic comparison
PCA, discriminant analysis, multidimensional scaling Multivariate methods
Worked R for every example in The Analysis of Biological Data Whitlock & Schluter examples

Filling the gaps.

Before you write it up

  • Analysis plan was written before data collection, and is in the repository
  • Replicate unit stated in words, and the model’s degrees of freedom agree with it
  • Non-independence is in the model
  • Residual diagnostic run, and the code committed
  • Candidate model set is small, and each model answers a stated question
  • Any AIC comparison of differing fixed effects used REML = FALSE
  • Estimate and 95% interval reported alongside every test
  • Effect size and both \(R^2\) values reported
  • Exact \(F\), both degrees of freedom, and exact \(p\)
  • Seed set and recorded for every resampling step
  • Package versions recorded
  • Any PCA states scaling, retention rule, variance per axis, and the sign convention
  • Any DMD reports rank, delay embedding, and error on data it did not fit
  • Measurement uncertainty reported separately, per Error and uncertainty analysis

Spot something to improve? Anyone in the lab can edit this guide. Edit this page suggest a change how edits work

Last updated

← Back to Lab Guide overview