The Box-Cox Family and Its Likelihood

Why Transform a Response?

A response transformation changes the scale on which a linear regression model describes its conditional mean and errors. The Box-Cox family provides a continuous set of power transformations, so the analysis can compare possibilities such as a reciprocal, a square root, no transformation, or a logarithm within one parameterized family. Its likelihood calculation is useful, but requires care: the conventional Gaussian calculation is generally a working or approximate likelihood, not an exact density for every nonzero power.

Let $Y_1,\ldots,Y_n$ be observed responses and let $\mathbf{x}_i$ be the predictor vector for observation $i$, including an intercept if appropriate. The Box-Cox family requires $Y_i>0$ for every observation. A zero or negative response cannot be inserted into the standard definition without changing the data or the model.

Definition: For a positive response $y$, the Box-Cox transformation with power $\lambda$ is $$ T_\lambda(y)= \begin{cases} \dfrac{y^\lambda-1}{\lambda}, & \lambda\ne 0,\\[6pt] \log y, & \lambda=0. \end{cases} $$ The value at $\lambda=0$ is defined by continuity. For $\lambda\ne 0$, the transformation is strictly increasing in $y>0$, with derivative $dT_\lambda(y)/dy=y^{\lambda-1}>0$.

Subtracting 1 and dividing by $\lambda$ makes this version of the power family approach the logarithm smoothly as $\lambda$ approaches zero. Some presentations use a geometric-mean-scaled variant to reduce dependence on the units used for $Y$. The transformation convention must therefore be stated when reporting a power analysis. The calculations below use the definition in the box.

The Zero-Power Limit and Transformation Values

For $y>0$, write $y^\lambda=e^{\lambda\log y}$. The expansion $e^a=1+a+o(a)$ as $a\to0$ gives

\lim_{\lambda\to0}\frac{y^\lambda-1}{\lambda} = \lim_{\lambda\to0}\frac{e^{\lambda\log y}-1}{\lambda} = \log y.

Thus, the logarithm is not an unrelated extra option: it is the continuous zero-power member of the family. A few other members are also familiar. $T_1(y)=y-1$, which differs from the untransformed response only by a constant shift; $T_{1/2}(y)=2(\sqrt y-1)$; and $T_{-1}(y)=1-1/y$.

Worked Example: Applying Four Box-Cox Powers

A hypothetical analyst records positive claim-processing costs, in a fixed monetary unit, and considers how one cost $y=4$ is transformed under $\lambda=1,\frac12,0,-1$. Calculate each transformed value.

Method. Substitute the same positive observation into the stated Box-Cox definition for each candidate power.

For $\lambda=1$, $$ T_1(4)=\frac{4^1-1}{1}=3. $$

For $\lambda=\frac12$, $$ T_{1/2}(4)=\frac{4^{1/2}-1}{1/2} =\frac{2-1}{1/2}=2. $$

For $\lambda=0$, use the defined limit: $$ T_0(4)=\log 4\approx1.3863. $$ For $\lambda=-1$, $$ T_{-1}(4)=\frac{4^{-1}-1}{-1} =\frac{1/4-1}{-1}=0.75. $$

The results are $3,2,1.3863,$ and $0.75$, respectively. The input is positive, as required. The common mistake this avoids is substituting $\lambda=0$ into the fraction, which would divide by zero; at zero the transformation is $\log y$.

The Gaussian Working Likelihood

For a candidate $\lambda$, the conventional Box-Cox regression calculation treats the transformed responses as approximately following a linear model:

T_\lambda(Y_i)=\mathbf{x}_i^{\mathsf T}\boldsymbol{\beta}+\varepsilon_i, \qquad \varepsilon_i\ \text{independent with}\ \varepsilon_i\sim N(0,\sigma^2).

Here the predictors are treated as fixed, the regression design has full column rank, and $\sigma^2>0$. This is a Gaussian model on the transformed scale used as a working approximation in the usual Box-Cox procedure. A crucial qualification is that, when $\lambda\ne0$, $T_\lambda(Y_i)$ cannot have an unrestricted normal distribution exactly: its range is restricted. The likelihood derived below is therefore not, in general, a normalized exact density for the original responses.

To see the conventional calculation, let $z_i=T_\lambda(y_i)$. The inverse transformation exists for positive $y_i$, and the absolute Jacobian derivative is $y_i^{\lambda-1}$. Applying the change-of-variables formula to the Gaussian working density gives the likelihood expression

L_{\mathrm{work}}(\boldsymbol{\beta},\sigma^2,\lambda) = (2\pi\sigma^2)^{-n/2} \exp\!\left\{-\frac{1}{2\sigma^2} \sum_{i=1}^n \big[T_\lambda(y_i)-\mathbf{x}_i^{\mathsf T}\boldsymbol{\beta}\big]^2\right\} \prod_{i=1}^n y_i^{\lambda-1}.

The product is the Jacobian contribution. In log form, it adds $(\lambda-1)\sum_i\log y_i$. This formula is commonly used to compare powers, but for nonzero $\lambda$ its Gaussian component assigns probability outside the transformation's range. Calling it a working likelihood keeps that limitation explicit.

Profiling Out the Regression Parameters

For each fixed $\lambda$, define $\operatorname{RSS}_\lambda$ as the residual sum of squares from the ordinary least-squares regression of $T_\lambda(y_i)$ on the design matrix:

\operatorname{RSS}_\lambda = \min_{\boldsymbol{\beta}} \sum_{i=1}^n \big[T_\lambda(y_i)-\mathbf{x}_i^{\mathsf T}\boldsymbol{\beta}\big]^2.

The minimizer is the least-squares estimate $\widehat{\boldsymbol{\beta}}_\lambda$. Substituting it into the likelihood leaves a function of $\sigma^2$ and $\lambda$. For fixed $\lambda$, differentiating the log-likelihood with respect to $\sigma^2$ gives the maximum-likelihood estimate $\widehat{\sigma}^2_\lambda=\operatorname{RSS}_\lambda/n$. This is the maximum-likelihood divisor $n$, not the unbiased-estimation divisor $n-p$, where $p$ is the number of fitted regression coefficients.

Choose a candidate power.
Transform every positive response using $T_\lambda$, using $\log y_i$ when $\lambda=0$.

Fit the transformed-response regression.
Use least squares to find $\widehat{\boldsymbol{\beta}}_\lambda$ and calculate $\operatorname{RSS}_\lambda$.

Profile out the error variance.
Set $\widehat{\sigma}^2_\lambda=\operatorname{RSS}_\lambda/n$ in the Gaussian working log-likelihood.

Compare candidate powers.
Evaluate the profile log-likelihood, including the Jacobian term, across candidate values of $\lambda$.

The resulting profile log-likelihood is

\ell_{\mathrm{work}}(\lambda) = -\frac{n}{2}\left[ \log(2\pi)+1+\log\!\left(\frac{\operatorname{RSS}_\lambda}{n}\right) \right] +(\lambda-1)\sum_{i=1}^n\log y_i.

The formula follows by taking logarithms of $L_{\mathrm{work}}$, replacing the residual sum of squares by its least-squares minimum, and then substituting $\widehat{\sigma}^2_\lambda=\operatorname{RSS}_\lambda/n$. For a fixed data set, the terms $-\frac n2[\log(2\pi)+1]$ are the same for all powers. They can be omitted when ranking candidate powers, but the RSS term and the Jacobian term cannot.

Worked Example: Comparing Three Profile Likelihood Values

For a hypothetical data set with three positive responses $1,2,4$, fit an intercept-only transformed-response model and compare the working profile log-likelihood at $\lambda=1,0,-1$. Report the values after dropping the common constant $-\frac{n}{2}[\log(2\pi)+1]$.

Method. For each power, transform the three observations, calculate the mean and residual sum of squares, and evaluate $$ q(\lambda)=-\frac{n}{2}\log\!\left(\frac{\operatorname{RSS}_\lambda}{n}\right) +(\lambda-1)\sum_i\log y_i. $$ Here $n=3$ and $\sum_i\log y_i=\log(1\cdot2\cdot4)=\log 8\approx2.0794$.

At $\lambda=1$, the transformed responses are $0,1,3$. Their mean is $4/3$, so $$ \operatorname{RSS}_1=(0-4/3)^2+(1-4/3)^2+(3-4/3)^2 =\frac{16+1+25}{9}=\frac{14}{3}. $$ Thus $\operatorname{RSS}_1/n=14/9$, and $$ q(1)=-\frac32\log(14/9)+0\cdot\log8\approx-0.6627. $$

At $\lambda=0$, the transformed values are $0,\log2,\log4$, with mean $\log2$. Hence $$ \operatorname{RSS}_0=(\log2)^2+0^2+(\log2)^2 =2(\log2)^2\approx0.9609. $$ Therefore $\operatorname{RSS}_0/3\approx0.3203$, and $$ q(0)=-\frac32\log(0.3203)-\log8\approx1.7077-2.0794=-0.3717. $$

At $\lambda=-1$, the transformed values are $0,\frac12,\frac34$, with mean $5/12$. Thus $$ \operatorname{RSS}_{-1} =\left(-\frac5{12}\right)^2+\left(\frac1{12}\right)^2+\left(\frac4{12}\right)^2 =\frac{42}{144}=\frac7{24}. $$ Now $\operatorname{RSS}_{-1}/3=7/72$, so $$ q(-1)=-\frac32\log(7/72)-2\log8\approx3.4961-4.1589=-0.6628. $$

Among these three candidates, $\lambda=0$ has the largest working profile value. The arithmetic is checked by using the intercept-only mean for each transformed response and retaining the Jacobian term. A common error is to compare only the residual sums of squares; the profile likelihood also includes $(\lambda-1)\sum_i\log y_i$.

What the Support Restriction Changes

For nonzero powers, the transformation does not range over all real numbers:

T_\lambda((0,\infty))= \begin{cases} (-1/\lambda,\infty), & \lambda>0,\\ (-\infty,-1/\lambda), & \lambda<0. \end{cases}

For example, when $\lambda>0$, $y^\lambda>0$ implies $T_\lambda(y)>-1/\lambda$. A normal random variable has positive probability below that bound. When $\lambda<0$, the transformed response is instead bounded above by $-1/\lambda$, and a normal random variable has positive probability above that bound. Thus an unrestricted normal distribution for $T_\lambda(Y_i)$ cannot be the exact conditional distribution of the transformed positive response for nonzero $\lambda$.

An exact model can instead specify a normal distribution truncated to the appropriate range. For example, write $\mu_i=\mathbf{x}_i^{\mathsf T}\boldsymbol{\beta}$, and let $\phi$ and $\Phi$ denote the standard normal density and distribution function. For $\lambda>0$, the exact density of $Y_i$ under a normal distribution truncated below at $-1/\lambda$ is

f_{Y_i}(y_i\mid\mathbf{x}_i) = \frac{ \phi\!\left((T_\lambda(y_i)-\mu_i)/\sigma\right) }{ \sigma\left[1-\Phi\!\left((-1/\lambda-\mu_i)/\sigma\right)\right] } y_i^{\lambda-1},\qquad y_i>0.

For $\lambda<0$, the truncation probability in the denominator is instead $\Phi((-1/\lambda-\mu_i)/\sigma)$. These normalizing probabilities depend on the model parameters, so they alter the likelihood and its profile. At $\lambda=0$, the transformation range is all real numbers; a normal model for $\log Y_i$ gives an exact lognormal model for $Y_i$, with no truncation factor.

Worked Example: Normalizing a Truncated Model

Suppose a hypothetical positive response follows an exact model with $\lambda=1$, and $T_1(Y)$, conditional on its predictors, has a $N(0,1)$ distribution truncated to its allowed range. Find the density at $y=1.5$, and compare it with the untruncated Gaussian working expression.

Method. Since $T_1(y)=y-1$, the transformed value is $T_1(1.5)=0.5$, and the support is $(-1,\infty)$. The probability that an untruncated $N(0,1)$ variable falls in this range is $$ P(Z>-1)=1-\Phi(-1)=\Phi(1)\approx0.8413. $$ The Jacobian is $y^{1-1}=1$. Consequently, the exact density is $$ f_Y(1.5)=\frac{\phi(0.5)}{\Phi(1)} =\frac{0.3521}{0.8413}\approx0.4184. $$

The untruncated working expression would give only $\phi(0.5)\approx0.3521$. It is not a normalized density on the allowed transformed range: the untruncated normal assigns only $0.8413$ probability to that range. Dividing by $0.8413$ corrects the normalization. The example shows why the ordinary Box-Cox profile is described as working or approximate, rather than as an exact likelihood under an unrestricted transformed normal model.

Interpreting a Selected Power

A selected power identifies a response scale that performs well under the specified regression and likelihood procedure. Common reference values give useful descriptions: a power near 1 corresponds to an untransformed response up to a constant shift, a power near 0 to a logarithmic response, a power near $1/2$ to a square-root-type response, and a power near $-1$ to a reciprocal-type response. These descriptions concern the scale of the modeled response, not a guarantee that the underlying data-generating process has exactly that form.

The regression coefficients describe conditional means on the transformed scale. In particular, a slope in a log-response model describes an additive change in the conditional mean of $\log Y$; exponentiating that change gives a multiplicative factor for the geometric scale. It does not, by itself, equal the change in the arithmetic conditional mean of $Y$. For other powers, interpretation likewise requires respecting the transformed scale and the inverse transformation. Selecting a power does not automatically establish a causal effect or make every regression assumption true.

The Box-Cox procedure is aligned with the SOA Exam SRM (Statistics for Risk Modeling) syllabus. In applied work, the response units, positivity requirement, chosen transformation convention, and working-likelihood qualification should be reported. A power should not be interpreted in isolation from the fitted model and the data.

Key takeaway: The Box-Cox profile combines transformed-scale fit, measured by $\operatorname{RSS}_\lambda$, with the change-of-variables Jacobian. Its conventional Gaussian form is a working likelihood for nonzero powers because the transformed response has restricted support; an exact truncated model has additional parameter-dependent normalizing factors.

Check Your Understanding

For each problem, use the stated Box-Cox convention and distinguish the working likelihood from an exact truncated model where relevant.

For a positive response $y=9$, calculate $T_{1/2}(y)$ and state why $y=0$ cannot be substituted in the standard family. (Answer: $T_{1/2}(9)=4$; the standard family requires $y>0$.)
For $\lambda=-1/2$, determine the range of $T_\lambda(Y)$ when $Y>0$. (Answer: $(-\infty,2)$.)
An intercept-only regression has $n=4$, $\operatorname{RSS}_{1/2}=8$, and $\sum_i\log y_i=\log16$. Calculate the profile quantity $q(\lambda)=-\frac n2\log(\operatorname{RSS}_\lambda/n)+(\lambda-1)\sum_i\log y_i$ at $\lambda=1/2$. (Answer: $-2\log2-\frac12\log16=-4\log2\approx-2.7726$.)
For $\lambda=1$, a model has $T_1(Y)\sim N(0,1)$ truncated to its allowed range. State the range and the probability used to normalize the density. (Answer: $(-1,\infty)$; the normalizing probability is $\Phi(1)\approx0.8413$.)
In a log-response regression, a one-unit predictor increase changes the fitted conditional mean of $\log Y$ by $0.2$. Give the corresponding multiplicative factor on the geometric scale and its percentage increase. (Answer: $e^{0.2}\approx1.2214$, or about a $22.14\%$ increase on the geometric scale.)

What You'll Learn