Yeo-Johnson Likelihood Inference in R

When the data are not normal enough.

Many statistical tools, from the t-test to linear models, assume roughly normal data. Real measurements are often skewed. A power transformation can make them more symmetric, but the classic Box-Cox family only works on strictly positive values.

Yeo and Johnson (Biometrika, 2000) proposed a family that accepts negative values, zero and positive values alike, which makes it the right tool for signed measurements such as the difference between two paired plants. I implemented it directly in R, with its logarithmic special cases:

One parameter, many shapes.

The second script studies the family for λ = 0, 0.5, 1, 1.5 and 2 and compares it with Box-Cox. At λ = 1 the transformation is the identity. Below 1 it compresses large positive values and stretches negative ones, which corrects right skew; above 1 it does the opposite and corrects left skew.

Computed from the transformation's formula, in the style of the original paper's figure. Curves are clipped at ±6.

Testing the limits on purpose.

The first script builds a deliberately asymmetric distribution: 30% standard normal and 70% Gamma(4, 1) shifted by −2, discretized on a 3,000-point grid, from which it draws 10,000 observations. It then picks the transformation parameter by bounded numerical optimization (L-BFGS-B) and compares three curves: the original mixture, a kernel density of the transformed sample, and the Gaussian fitted to it.

The point is to see how far a single power transformation goes on a shape it was not designed for, and where it stops helping. That is what motivates the more flexible Yeo-Johnson family in the next step.

Back to 1876.

Darwin grew 15 pairs of plants, one cross-fertilized and one self-fertilized in each pair, and measured their heights. The question is simple: do cross-fertilized plants grow taller? The data are the 15 differences, two of them negative, the same ones the demo at the top transforms.

Under a normal model for the transformed values, the log-likelihood includes the Jacobian of the transformation, the term that connects the model to the original measurements:

The script optimizes λ first, then refines λ, μ and σ² jointly, and computes their covariance matrix from the numerical Hessian. The estimate is λ̂ = 1.305, the same value Yeo and Johnson report for this data. The transformation fixes the shape:

After the transformation the data are almost perfectly symmetric, normality is far more plausible, and the evidence that cross-fertilized plants grow taller gets stronger.

yeo_johnson <- function(x, lambda) {  y <- numeric(length(x))  pos_values <- x >= 0  neg_values <- x < 0  # negative values: a modified power transformation, continuous across zero  y[neg_values] <- if (lambda == 2) {    -log(-x[neg_values] + 1)  } else {    -((-x[neg_values] + 1)^(2 - lambda) - 1) / (2 - lambda)  }  # positive values: the Box-Cox form, shifted by one  y[pos_values] <- if (lambda == 0) {    log(x[pos_values] + 1)  } else {    ((x[pos_values] + 1)^lambda - 1) / lambda  }  return(y)}log_likelihood <- function(lambda, mu, sigma2, x) {  transformed_x <- yeo_johnson(x, lambda)  n <- length(x)  # Gaussian likelihood of the transformed data plus the Jacobian of the transformation  log_likelihood_value <- -(n / 2) * log(2 * pi) - (n / 2) * log(sigma2) -    (1 / (2 * sigma2)) * sum((transformed_x - mu)^2) +    (lambda - 1) * sum(sign(x) * log(abs(x) + 1))  return(-log_likelihood_value)   # negative, for minimization with optim()}

Is the transformation worth it?

A transformation is a modeling choice, so it has to earn its place. A likelihood-ratio test compares the optimized model with λ = 1, the untransformed data: the statistic is 3.92 against a χ² distribution with one degree of freedom, a p-value of about 0.048. The improvement is real, at the edge of the usual 5% threshold.

The script closes by returning to the original scale: it back-transforms a lower 0.01 quantile for the mean difference, based on a Student-t reference, so that the conclusion is stated in the units Darwin measured.