{"word_count": 435, "character_count": 3732}
A few days ago, a former student of mine, David, came back to me about Box-Cox tests in linear models. It made me look more carefully at the test, and I do not understand what is computed, to be honest. Let us start with something simple, like a linear simple regression, i.e.Let us introduced - as suggested in Box &#38; Cox (1964) - the following family of (power) transformationson the variable of interest. Then assume thatAs mentioned in Chapter 14 of Davidson &#38; MacKinnon (1993) - in French - the log-likelihood of this model (assuming that observations are independent, with distribution ) can be writtenWe can then use profile-likelihood techniques (see here) to derive the optimal transformation.This can be done in R extremely simply,&#62; library(MASS)&#62; boxcox(lm(dist~speed,data=cars),lambda=seq(0,1,by=.1))we then get the following graph,If we look at the code of the function, it is based on the QR decomposition of the matrix (since we assume that  is a full-rank matrix). More precisely,  where  is a  matrix,  is a  orthonornal matrix, and is a  upper triangle matrix. It might be convenient to use this matrix since, for instance, .  Thus, we do have an upper triangle system of equations.&#62; X=lm(dist~speed,data=cars)$qrThe code used to get the previous graph is (more or less) the following,&#62; g=function(x,lambda){+ y=NA+ if(lambda!=0){y=(x^lambda-1)/lambda}+ if(lambda==0){y=log(x)}+ return(y)} &#62; n=nrow(cars)&#62; X=lm(dist~speed,data=cars)$qr&#62; Y=cars$dist&#62; logv=function(lambda){+ -n/2*log(sum(qr.resid(X, g(Y,lambda)/+ exp(mean(log(Y)))^(lambda-1))^2))}&#62; L=seq(0,1,by=.05)&#62; LV=Vectorize(logv)(L)&#62; points(L,LV,pch=19,cex=.85,col="red")As we can see (with those red dots) we can reproduce the R graph. But it might not be consistent with other techniques (and functions described above). For instance, we can plot the profile likelihood function, &#62; logv=function(lambda){+ s=summary(lm(g(dist,lambda)~speed,+ data=cars))$sigma+ e=lm(g(dist,lambda)~speed,data=cars)$residuals+ -n/2*log(2 * pi)-n*log(s)-.5/s^2*(sum(e^2))++ (lambda-1)*sum(log(Y))+ }&#62; L=seq(0,1,by=.01)&#62; LV=Vectorize(logv)(L)&#62; plot(L,LV,type="l",ylab="")&#62; (maxf=optimize(logv,0:1,maximum=TRUE))$maximum[1] 0.430591$objective[1] -197.6966&#62; abline(v=maxf$maximum,lty=2)The good point is that the optimal value of  is the same as the one we got before. The only problem is that the -axis has a different scale. And using profile likelihood techniques to derive a confidence interval will give us different results (with a larger confidence interval than the one given by the standard function),&#62; ic=maxf$objective-qchisq(.95,1)&#62; #install.packages("rootSolve")&#62; library(rootSolve)&#62; f=function(x)(logv(x)-ic)&#62; (lower=uniroot(f, c(0,maxf$maximum))$root)[1] 0.1383507&#62; (upper=uniroot(f, c(maxf$maximum,1))$root)[1] 0.780573&#62; segments(lower,ic,upper,ic,lwd=2,col="red")Actually, it possible to rewrite the log-likelihood as(let us just get rid of the constant), whereHere, it becomes&#62; logv=function(lambda){+ e=lm(g(dist,lambda)~speed,data=cars)$residuals+ elY=(exp(mean(log(Y))))+ -n/2*log(sum((e/elY^lambda)^2))+ }&#62;&#62; L=seq(0,1,by=.01)&#62; LV=Vectorize(logv)(L)&#62; plot(L,LV,type="l",ylab="")&#62; optimize(logv,0:1,maximum=TRUE)$maximum[1] 0.430591$objective[1] -47.73436with again the same optimal value for , and the same confidence interval, since the function is the same, up to some additive constant.So we have been able to derive the optimal transformation according to Box-Cox transformation, but so far, the confidence interval is not the same (it might come from the fact that here we substituted an estimator to the unknown parameter .