6.3 The Box–Cox method
Rather than guess a transformation, we can let the data choose one from the family of power transformations \(y^\lambda \). Box and Cox (1964) made this precise. Some care is needed, because the residual sums of squares from fits to \(y^\lambda \) for different \(\lambda \) are on different scales and cannot be compared directly. The remedy is to rescale.
Definition 6.3.1. For \(y_i > 0\), let \(\dot y = \bigl (\prod _{i=1}^n y_i\bigr )^{1/n}\) be the geometric mean of the responses. The Box–Cox transformation is \[ y^{(\lambda )} = \begin {cases} \dfrac {y^\lambda - 1}{\lambda \,\dot y^{\,\lambda - 1}} , & \lambda \neq 0, \\[3mm] \dot y\ln y , & \lambda = 0 . \end {cases} \]
Subtracting \(1\) and dividing by \(\lambda \) does not change the fit, apart from scaling; it makes \(y^{(\lambda )}\) vary continuously with \(\lambda \), with the \(\lambda = 0\) case as the limit (Problem 6.5.2). Dividing by \(\dot y^{\lambda - 1}\) puts every \(\lambda \) on the scale of \(y\) itself, so that the residual sums of squares can be compared.
Theorem 6.3.2. Suppose that for some \(\lambda \), \(y_i^{(\lambda )} = \mathbf {x}_i^\top \boldsymbol \beta + \varepsilon _i\) with \(\varepsilon _i\) independent \(N(0, \sigma ^2)\). Let \(SS_{\mathrm {Res}}(\lambda )\) be the residual sum of squares from the least-squares fit of \(y^{(\lambda )}\) on \(\mathbf {X}\). Then, up to an additive constant, the log-likelihood of the observed \(y_i\), maximised over \(\boldsymbol \beta \) and \(\sigma ^2\) for fixed \(\lambda \), is \[ L(\lambda ) = -\frac {n}{2}\ln SS_{\mathrm {Res}}(\lambda ) . \] So the maximum-likelihood estimate \(\hat \lambda \) is the value that minimises \(SS_{\mathrm {Res}}(\lambda )\).
Proof. The density of \(w_i = y_i^{(\lambda )}\) is normal; the density of \(y_i\) is obtained by the change of variables formula, multiplying by the Jacobian \(|dw_i/dy_i|\). For \(\lambda \neq 0\), \[ \frac {dw_i}{dy_i} = \frac {\lambda y_i^{\lambda - 1}}{\lambda \,\dot y^{\,\lambda - 1}} = \left (\frac {y_i}{\dot y}\right )^{\lambda - 1} , \] and for \(\lambda = 0\) it is \(\dot y/y_i\), the same formula. The product over \(i\) is \[ \prod _{i=1}^n\left (\frac {y_i}{\dot y}\right )^{\lambda - 1} = \left (\frac {\prod y_i}{\dot y^{\,n}}\right )^{\lambda - 1} = 1^{\lambda - 1} = 1 , \] by the definition of \(\dot y\). So the likelihood of the \(y_i\) equals the normal likelihood of the \(w_i\): \[ \ln \ell = -\frac {n}{2}\ln (2\pi \sigma ^2) - \frac {1}{2\sigma ^2}\sum _{i=1}^n\bigl (w_i - \mathbf {x}_i^\top \boldsymbol \beta \bigr )^2 . \] For fixed \(\sigma ^2\) this is maximised over \(\boldsymbol \beta \) by least squares, leaving \(-\frac {n}{2}\ln (2\pi \sigma ^2) - SS_{\mathrm {Res}}(\lambda )/(2\sigma ^2)\). Exactly as in Theorem 2.6.1, this is maximised at \(\hat \sigma ^2 = SS_{\mathrm {Res}}(\lambda )/n\), where its value is \[ -\frac {n}{2}\ln SS_{\mathrm {Res}}(\lambda ) + \frac {n}{2}\ln n - \frac {n}{2}\ln (2\pi ) - \frac {n}{2} . \] The last three terms do not depend on \(\lambda \). □
Without the factor \(\dot y^{\lambda - 1}\) the Jacobian would be \(\prod y_i^{\lambda - 1}\), which depends on \(\lambda \), and minimising the raw residual sum of squares would give the wrong answer. The normalisation absorbs the Jacobian.
A confidence interval for \(\lambda \)
By the large-sample theory of likelihood ratios, \(2\bigl [L(\hat \lambda ) - L(\lambda )\bigr ]\) is approximately \(\chi ^2_1\) when \(\lambda \) is the true value. So an approximate \(100(1 - \alpha )\%\) confidence interval is the set of \(\lambda \) with \(L(\lambda ) \ge L(\hat \lambda ) - \frac 12\chi ^2_{\alpha ,1}\), that is, with \[ SS_{\mathrm {Res}}(\lambda ) \le SS^* = SS_{\mathrm {Res}}(\hat \lambda )\, e^{\chi ^2_{\alpha ,1}/n} . \] In practice \(SS_{\mathrm {Res}}(\lambda )\) is computed for a grid of values, say from \(-2\) to \(2\), and plotted. One then chooses a simple \(\lambda \) inside the interval, such as \(0\), \(\frac 12\) or \(1\), in preference to the exact minimiser: a transformation such as \(y^{0.43}\) is hard to interpret, and the data rarely distinguish it from \(\sqrt {y}\). If \(1\) is inside the interval, there is no evidence that any transformation is needed.
Example 6.3.3. Apply the Box–Cox method to the bacteria data of Example 6.2.2.
Solution. The geometric mean is \(\dot y = \exp (\overline {\ln y}) = e^{3.72527} = 41.482\). For each \(\lambda \), form \(y^{(\lambda )}\) and regress it on \(t\). For instance, at \(\lambda = 0.5\), \(y^{(0.5)} = (\sqrt {y} - 1)/(0.5 \times 41.482^{-0.5}) = 12.881(\sqrt {y} - 1)\), and at \(\lambda = 0\), \(y^{(0)} = 41.482\ln y\), whose residual sum of squares is \(41.482^2 \times 0.03069 = 52.81\). At \(\lambda = 1\), \(y^{(1)} = y - 1\), which has the same residuals as \(y\). The results:
| \(\lambda \) | \(-1\) | \(-0.5\) | \(-0.25\) | \(0\) | \(0.25\) | \(0.5\) | \(1\) |
| \(SS_{\mathrm {Res}}(\lambda )\) | 3223.4 | 609.0 | 172.1 | 52.8 | 194.5 | 620.4 | 2871.6 |
A finer search puts the minimum at \(\hat \lambda = -0.01\), with \(SS_{\mathrm {Res}}(\hat \lambda ) = 52.44\). With \(n = 8\) and \(\chi ^2_{0.05,1} = 3.841\), \[ SS^* = 52.44\, e^{3.841/8} = 52.44 \times 1.616 = 84.76 , \] and \(SS_{\mathrm {Res}}(\lambda ) \le 84.76\) for \(-0.14 \le \lambda \le 0.11\). The interval contains \(0\) and excludes both \(\frac 12\) and \(1\): the log transformation is the right choice, and neither the square root nor the untransformed response is acceptable. This confirms the choice made in Example 6.2.2 on physical grounds.
Remark (Transforming the regressors). Regressors can be transformed as well (Example 6.2.1). A transformation of \(x\) changes the shape of the fitted relationship but not the error structure, so it is the remedy for curvature when the variance is already constant, and a transformation of \(y\) is the remedy when it is not. The partial residual plots of Section 4 are the main guide to the form in which each regressor should enter.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.