Usually transformation to numeric data D = { X i , Y i } i = 1 n \mathcal{D}=\{X_i,Y_i\}_{i=1}^n D = { X i , Y i } i = 1 n is necessary in regression analysis, usually to stablize variance . Box-Cox is the most important method.
Y ∗ = Y λ − 1 λ e . g . = { Y ∗ ∼ Y , λ = 1 Y ∗ ∼ Y , λ = 0.5 Y ∗ ∼ ln Y , λ = 0 Y ∗ ∼ 1 / Y , λ = − 1 \begin{align}
Y^*=\dfrac{Y^\lambda -1}{\lambda }\mathrm{e.g.}=\begin{cases}
Y^*\sim Y &,\lambda =1\\
Y^*\sim \sqrt{Y} &,\lambda =0.5\\
Y^*\sim \ln Y &,\lambda =0\\
Y^*\sim 1\big/ Y &,\lambda =-1
\end{cases}
\end{align} Y ∗ = λ Y λ − 1 e.g. = ⎩ ⎨ ⎧ Y ∗ ∼ Y Y ∗ ∼ Y Y ∗ ∼ ln Y Y ∗ ∼ 1 / Y , λ = 1 , λ = 0.5 , λ = 0 , λ = − 1
with linear model
Y ∗ = X ′ β ∗ + ε ∗ , ε ∗ ∼ N ( 0 , σ 2 ) \begin{align}
Y^*=X'\beta^* +\varepsilon^* ,\quad \varepsilon^* \sim N(0,\sigma ^2)
\end{align} Y ∗ = X ′ β ∗ + ε ∗ , ε ∗ ∼ N ( 0 , σ 2 )
where X = ( 1 , X 1 , … , X p ) ′ X=(1,X_1,\ldots,X_p)' X = ( 1 , X 1 , … , X p ) ′ , β ∗ = ( β 0 , β 1 , … , β p ) ′ \beta ^*=(\beta _0,\beta _1,\ldots,\beta _p)' β ∗ = ( β 0 , β 1 , … , β p ) ′
Likelihood function expressed in D = { X i , Y i ∗ } i = 1 n \mathcal{D}=\{X_i,Y_i^*\}_{i=1}^n D = { X i , Y i ∗ } i = 1 n :
L ( β ∗ , σ 2 ; λ ) = 1 ( 2 π σ 2 ) n / 2 exp [ − 1 2 σ 2 ∑ i = 1 n ( Y i ∗ − X i ′ β ∗ ) 2 ] ∣ J ( ∂ Y ∗ ∂ Y ) ∣ \begin{align}
L(\beta^* ,\sigma ^2;\lambda )=\dfrac{1}{(2\pi\sigma ^2)^{n/2}}\exp\left[ -\dfrac{1}{2\sigma ^2}\sum_{i=1}^n\left( Y_i^*-X_i'\beta^* \right)^2 \right] \left|J\left(\dfrac{\partial^{} Y^*}{\partial Y^{}}\right)\right|
\end{align} L ( β ∗ , σ 2 ; λ ) = ( 2 π σ 2 ) n /2 1 exp [ − 2 σ 2 1 i = 1 ∑ n ( Y i ∗ − X i ′ β ∗ ) 2 ] J ( ∂ Y ∂ Y ∗ )
where the Jacobi matrix could be denoted in Geometric Mean G M ( Y ) = ∏ i = 1 n Y i 1 / n \mathrm{GM}(Y )=\prod_{i=1}^n Y_i^{1/n} GM ( Y ) = ∏ i = 1 n Y i 1/ n
∣ J ( ∂ Y ∗ ∂ Y ) ∣ = ∏ i = 1 n Y i λ − 1 G M ( Y ) n ( λ − 1 ) \begin{align}
\left|J\left(\dfrac{\partial^{} Y^*}{\partial Y^{}}\right)\right|=\prod_{i=1}^nY_i^{\lambda -1}\mathrm{GM}\left( Y \right) ^{n(\lambda -1)}
\end{align} J ( ∂ Y ∂ Y ∗ ) = i = 1 ∏ n Y i λ − 1 GM ( Y ) n ( λ − 1 )
MLE estimators are similar:
β ^ ∗ = ( X ′ X ) − 1 X ′ Y ∗ σ ^ n 2 = 1 n ∑ i = 1 n ( Y i ∗ − Y ˉ ∗ ) \begin{align}
\hat{\beta }^*=&(X'X)^{-1}X'Y^*\\
\hat{\sigma }^2_n=&\dfrac{1}{n}\sum_{i=1}^n(Y_i^*-\bar{Y}^*)
\end{align} β ^ ∗ = σ ^ n 2 = ( X ′ X ) − 1 X ′ Y ∗ n 1 i = 1 ∑ n ( Y i ∗ − Y ˉ ∗ )
Subtitute MLE estimators back to (log-)likelihood:
log L ( β , σ 2 ; λ ) = ℓ ( λ ) = − n 2 log σ ^ n 2 G M ( Y ) 2 ( λ − 1 ) + c o n s t \begin{align}
\log L(\beta ,\sigma ^2;\lambda )=\ell(\lambda )=-\dfrac{n}{2}\log \dfrac{\hat{\sigma }^2_n}{\mathrm{GM}(Y)^{2(\lambda -1)} }+\mathrm{const}
\end{align} log L ( β , σ 2 ; λ ) = ℓ ( λ ) = − 2 n log GM ( Y ) 2 ( λ − 1 ) σ ^ n 2 + const
By plotting ℓ ( λ ) \ell(\lambda ) ℓ ( λ ) v.s. λ \lambda λ we could locate a λ \lambda λ both appropriate in interpretability and likelihood maximization. Here's the example code from r::MASS, which indicating a selection of λ = 0 \lambda =0 λ = 0 for logarithm transformation.
library(MASS)
boxcox(Volume ~ log(Height) + log(Girth), data = trees,
lambda = seq(-0.25, 0.25, length = 10))