Testing isometry and comparing slopes

R
allometry
body size
hypothesis testing
ecology tutorial
Is the exponent really 3/4? Testing a hypothesised allometric slope, comparing slopes between groups, and separating elevation from a shift along the line.
Author

Tidy Ecology

Published

2026-07-27

Fitting an exponent is rarely the end of an allometric analysis. The questions that follow are comparisons: is the slope different from the isometric value, is it different from 3/4, do the two sexes scale the same way, does the island population sit above the mainland one. Each of these is a hypothesis about a line, and each has a version that answers the wrong question convincingly.

This post covers three of them. The first is testing a slope against a hypothesised value, where the choice of estimator changes the null rejection rate by an order of magnitude. The second is testing whether several groups share a slope. The third is the distinction that gets lost most often: two groups can sit on the same line and still have very different mean trait values.

Testing a slope against a fixed value

Suppose the true exponent is 0.75 and you want to test that hypothesis. There is a neat exact test that comes straight out of the identity in the previous post. Build two derived variables from the data:

sma_slope <- function(x, y) sign(cor(x, y)) * sd(y) / sd(x)

pitman_test <- function(x, y, b0) {
  resid_score <- y - b0 * x
  fit_score   <- y + b0 * x
  cor.test(resid_score, fit_score)
}

The covariance of those two scores is var(y) - b0^2 * var(x), which is zero exactly when the population standardised major axis slope equals b0. So a plain correlation test on the residual and fitted scores is a test of the slope, and it is exact under bivariate normality rather than asymptotic.

Now the comparison that matters. Simulate data where the null is true, with body size measured imperfectly and the two error variances set so the standardised major axis condition holds, and count how often each approach rejects.

b_null <- 0.75
n_ind <- 60
n_mc <- 3000
sd_ex <- 0.30

set.seed(43811)
res <- t(vapply(seq_len(n_mc), function(i) {
  xt <- rnorm(n_ind, 3, 0.8)
  x <- xt + rnorm(n_ind, 0, sd_ex)
  y <- b_null * xt + rnorm(n_ind, 0, b_null * sd_ex)
  co <- summary(lm(y ~ x))$coefficients[2, ]
  p_ols <- 2 * pt(-abs((co[1] - b_null) / co[2]), n_ind - 2)
  c(p_ols, pitman_test(x, y, b_null)$p.value, co[1], sma_slope(x, y))
}, numeric(4)))
c(ols_reject = mean(res[, 1] < 0.05), pitman_reject = mean(res[, 2] < 0.05),
  mean_ols = mean(res[, 3]), mean_sma = mean(res[, 4]))
   ols_reject pitman_reject      mean_ols      mean_sma 
   0.48833333    0.04966667    0.65796523    0.75147938 

The null is true. The least squares t-test rejects it 48.8 percent of the time at a nominal five percent level, and the reason is in the last two numbers: the least squares slope averages 0.658 because measurement error in body size attenuates it, so the test is not testing a false hypothesis, it is testing a true one with a biased estimator. The Pitman test holds at 0.050.

Before concluding that the standardised major axis test is the safe default, sweep the amount of error in x.

sweep_x <- c(0, 0.10, 0.20, 0.30, 0.40)
sw <- do.call(rbind, lapply(seq_along(sweep_x), function(k) {
  su <- sweep_x[k]
  set.seed(9000 + k)
  o <- t(vapply(seq_len(1500), function(i) {
    xt <- rnorm(n_ind, 3, 0.8)
    x <- xt + rnorm(n_ind, 0, su)
    y <- b_null * xt + rnorm(n_ind, 0, sqrt((b_null * su)^2 + 0.12^2))
    co <- summary(lm(y ~ x))$coefficients[2, ]
    c(2 * pt(-abs((co[1] - b_null) / co[2]), n_ind - 2),
      pitman_test(x, y, b_null)$p.value)
  }, numeric(2)))
  data.frame(sd_error_x = su,
             ols = round(mean(o[, 1] < 0.05), 3),
             pitman_sma = round(mean(o[, 2] < 0.05), 3))
}))
sw
  sd_error_x   ols pitman_sma
1        0.0 0.045      0.132
2        0.1 0.071      0.091
3        0.2 0.205      0.054
4        0.3 0.441      0.051
5        0.4 0.664      0.049

Both columns should read 0.05 and neither does everywhere. Least squares climbs from 0.045 to 0.664 as body size gets noisier. The standardised major axis test goes the other way: at 0.132 when x is measured perfectly, settling near nominal only once the x error is large enough to balance the equation error in y.

That is the same condition as the previous post, showing up as a test size rather than as a bias. There is no estimator here that is safe in both regimes, and the sweep is the honest way to present the choice: you are picking which failure mode you are exposed to.

Line chart with error in x on the horizontal axis and rejection rate on the vertical. The OLS line rises steeply from near five percent. The SMA line starts high and falls towards five percent.
Figure 1: Rejection rate of a true null hypothesis about the slope, as measurement error in body size increases. Neither test holds its nominal size across the whole range, and they fail at opposite ends.

Do the groups share a slope?

The next question is comparative. Three populations, one relationship each, and you want to know whether a single exponent describes all of them. For the standardised major axis this is a likelihood ratio test, and it is short enough to write out.

At a candidate common slope b, each group has a correlation between its residual and fitted scores. That correlation is zero only at that group’s own standardised major axis slope, so forcing a common b makes them non-zero, and the total penalty is the statistic.

lr_common <- function(xs, ys) {
  pen <- function(b) {
    -sum(vapply(seq_along(xs), function(i) {
      rr <- cor(ys[[i]] - b * xs[[i]], ys[[i]] + b * xs[[i]])
      length(xs[[i]]) * log(1 - rr^2)
    }, numeric(1)))
  }
  own <- vapply(seq_along(xs), function(i) sma_slope(xs[[i]], ys[[i]]), numeric(1))
  op <- optimize(pen, c(min(own) * 0.5, max(own) * 1.5))
  list(stat = op$objective, common = op$minimum,
       p = pchisq(op$objective, length(xs) - 1, lower.tail = FALSE))
}

gen_group <- function(n, b, m, sd_size = 0.55, sd_err = 0.16) {
  xt <- rnorm(n, m, sd_size)
  list(x = xt + rnorm(n, 0, sd_err * 0.6),
       y = b * xt + rnorm(n, 0, b * sd_err * 0.6) + rnorm(n, 0, sd_err))
}

The chi-squared reference distribution is asymptotic, so check it rather than trusting it.

set.seed(4382)
null_p <- replicate(1500, {
  g <- lapply(1:3, function(k) gen_group(45, 0.75, 3))
  lr_common(lapply(g, `[[`, "x"), lapply(g, `[[`, "y"))$p
})
set.seed(4383)
alt_p <- replicate(1000, {
  g <- list(gen_group(45, 0.75, 3), gen_group(45, 0.75, 3), gen_group(45, 0.92, 3))
  lr_common(lapply(g, `[[`, "x"), lapply(g, `[[`, "y"))$p
})
set.seed(4384)
naive_p <- replicate(1000, {
  g <- lapply(1:3, function(k) gen_group(45, 0.75, 3))
  d <- do.call(rbind, lapply(seq_along(g), function(i)
    data.frame(x = g[[i]]$x, y = g[[i]]$y, grp = factor(i))))
  anova(lm(y ~ x * grp, d))["x:grp", "Pr(>F)"]
})
c(lr_null = mean(null_p < 0.05), lr_power = mean(alt_p < 0.05),
  naive_interaction_null = mean(naive_p < 0.05))
               lr_null               lr_power naive_interaction_null 
            0.05266667             0.60400000             0.05300000 

The likelihood ratio test sits at 0.053 under the null with three groups of 45, close enough to nominal to use, and detects a slope difference of 0.75 against 0.92 60 percent of the time. The interaction term from lm(y ~ x * grp) is also well calibrated, at 0.053.

So this is not a case where the familiar approach is broken. It is a case where the two approaches answer different questions. The interaction term compares least squares slopes, which are all attenuated by whatever measurement error is present; the likelihood ratio test compares standardised major axis slopes. If the groups differ in how precisely body size was measured, the attenuation differs between them, and the interaction term can report a slope difference that is entirely an artefact of field protocol. That is the reason to prefer the structural test, not a difference in power.

Elevation is not the same as position along the line

Suppose the common slope test passes and the groups share an exponent. There are still two ways they can differ, and conflating them is the most common error in this corner of the literature.

Construct two groups that lie on exactly the same line, differing only in the body sizes they contain.

set.seed(43822)
b_common <- 0.75
n_grp <- 70
sd_ex2 <- 0.1733
make_group <- function(mean_size) {
  xt <- rnorm(n_grp, mean_size, 0.45)
  list(x = xt + rnorm(n_grp, 0, sd_ex2),
       y = b_common * xt + rnorm(n_grp, 0, b_common * sd_ex2))
}
grp_a <- make_group(2.6)
grp_b <- make_group(3.4)

fit_common <- lr_common(list(grp_a$x, grp_b$x), list(grp_a$y, grp_b$y))
bc <- fit_common$common
res_score <- function(g) g$y - bc * g$x
fit_score <- function(g) g$y + bc * g$x

t_raw <- t.test(grp_a$y, grp_b$y)
t_elev <- t.test(res_score(grp_a), res_score(grp_b))
t_shift <- t.test(fit_score(grp_a), fit_score(grp_b))
c(common_slope = bc, p_raw_means = t_raw$p.value,
  p_elevation = t_elev$p.value, p_shift_along = t_shift$p.value)
 common_slope   p_raw_means   p_elevation p_shift_along 
 7.335990e-01  1.102082e-12  3.932045e-01  9.470700e-13 

The two groups were generated from one relationship with one exponent. The only thing that differs is mean body size.

Compare the raw trait means and you get a p-value of 1.1e-12, which reads as a large and highly significant group difference. Test elevation properly, by comparing residual scores about the common line, and the p-value is 0.393: no difference, which is correct, because there is none. Test position along the line and you recover 9.5e-13, which is the real and rather less exciting finding that one group contains larger animals.

Scatter with two coloured groups of points lying along the same straight line, one group occupying the lower left and the other the upper right, with a single fitted line through both.
Figure 2: Two groups on a shared allometric line. They differ in where they sit along the line, not in how far they sit above or below it, and a comparison of raw trait means cannot tell those apart.

The biological readings are not close. A shift in elevation says that at a given body size, one group has more of the trait: a grade difference, the kind of thing that suggests an adaptive shift. A shift along the axis says only that one group is bigger, and the trait follows along the same rule it always did. Report the second as the first and you have manufactured a result.

Where to go next

The sequence for a grouped allometric comparison runs in a fixed order, because each step depends on the one before. Test for a common slope first; if the slopes differ, elevation is not defined and there is nothing further to compare. If they share a slope, test elevation about the common line. Then, separately, test position along the line, and report it as what it is.

Honest limits

The Pitman test is exact under bivariate normality and only approximate otherwise, and the likelihood ratio test above leans on a chi-squared reference that was checked here for three groups of 45 and no smaller. The whole apparatus tests slopes of the fitted line, which is only the biological exponent under the error model the line assumes. And a common slope test that fails to reject is not evidence that the slopes are equal: with 45 animals per group the power against a difference of 0.17 was 60 percent, so a null result is compatible with a difference that would matter.

References

Warton, Weber 2002 Biometrical Journal 44(2):161-174 (10.1002/1521-4036(200203)44:2<161::AID-BIMJ161>3.0.CO;2-N)

Warton, Wright, Falster, Westoby 2006 Biological Reviews 81(2):259-291 (10.1017/S1464793106007007)

Newsletter

Get new tutorials by email

New R and QGIS tutorials for ecologists, straight to your inbox. No spam; unsubscribe anytime.

By subscribing you agree to receive these emails and confirm your address once. See the privacy policy.