Thanks John.
?boxcox says:
*************************
Arguments
object
a formula or fitted model object. Currently only lm and aov objects are handled.
*************************
I read that as saying that
boxcox(lm(z+1 ~ 1),...)
should run without error. But it didn't. And perhaps here's why:
BoxCoxLambda <- function(z){
b <- MASS:::boxcox.lm(lm(z+1 ~ 1), lambda = seq(-5, 5, length.out 61),
plotit = FALSE)
b$x[which.max(b$y)] # best lambda
}
> lambdas <- apply(dd,2 , BoxCoxLambda)
Error in NextMethod() : 'NextMethod' called from an anonymous function
and, indeed, ?UseMethod says:
"NextMethod should not be called except in methods called by UseMethod
or from internal generics (see InternalGenerics). In particular it
will not work inside anonymous calling functions (e.g.,
get("print.ts")(AirPassengers))."
****BUT ****
BoxCoxLambda <- function(z){
b <- MASS:::boxcox(z+1 ~ 1, lambda = seq(-5, 5, length.out = 61),
plotit = FALSE)
b$x[which.max(b$y)] # best lambda
}
> lambdas <- apply(dd,2 , BoxCoxLambda)
> lambdas
[1] 0.1666667 0.1666667
The identical lambdas do not seem right to me; nor do I understand why
boxcox.lm apparently throws the error while boxcox.formula does not
(it also calls NextMethod()) So I would welcome clarification to clear
my clogged (cerebral) sinuses. :-)
Best,
Bert
On Sat, Jul 8, 2023 at 11:25?AM John Fox <jfox at mcmaster.ca>
wrote:>
> Dear Ron and Bert,
>
> First (and without considering why one would want to do this, e.g.,
> adding a start of 1 to the data), the following works for me:
>
> ------ snip ------
>
> > library(MASS)
>
> > BoxCoxLambda <- function(z){
> + b <- boxcox(z + 1 ~ 1,
> + lambda = seq(-5, 5, length.out = 101),
> + plotit = FALSE)
> + b$x[which.max(b$y)]
> + }
>
> > mrow <- 500
> > mcol <- 2
> > set.seed(12345)
> > dd <- matrix(rgamma(mrow*mcol, shape = 2, scale = 5), nrow = mrow,
ncol > + mcol)
>
> > dd1 <- dd[, 1] # 1st column of dd
> > res <- boxcox(lm(dd1 + 1 ~ 1), lambda = seq(-5, 5, length.out =
101),
> plotit
> + = FALSE)
> > res$x[which.max(res$y)]
> [1] 0.2
>
> > apply(dd, 2, BoxCoxLambda, simplify = TRUE)
> [1] 0.2 0.2
>
> ------ snip ------
>
> One could also use the powerTransform() function in the car package,
> which in this context transforms towards *multi*normality:
>
> ------ snip ------
>
> > library(car)
> Loading required package: carData
>
> > powerTransform(dd + 1)
> Estimated transformation parameters
> Y1 Y2
> 0.1740200 0.2089925
>
> I hope this helps,
> John
>
> --
> John Fox, Professor Emeritus
> McMaster University
> Hamilton, Ontario, Canada
> web: https://www.john-fox.ca/
>
> On 2023-07-08 12:47 p.m., Bert Gunter wrote:
> > Caution: This email may have originated from outside the organization.
Please exercise additional caution with any links and attachments.
> >
> >
> > No, I'm afraid I'm wrong. Something went wrong with my R
session and gave
> > me incorrect answers. After restarting, I continued to get the same
error
> > as you did with my supposed "fix." So just ignore what I
said and sorry for
> > the noise.
> >
> > -- Bert
> >
> > On Sat, Jul 8, 2023 at 8:28?AM Bert Gunter <bgunter.4567 at
gmail.com> wrote:
> >
> >> Try this for your function:
> >>
> >> BoxCoxLambda <- function(z){
> >> y <- z
> >> b <- boxcox(y + 1 ~ 1,lambda = seq(-5, 5, length.out = 61),
plotit > >> FALSE)
> >> b$x[which.max(b$y)] # best lambda
> >> }
> >>
> >> ***I think*** (corrections and clarification strongly welcomed!)
that `~`
> >> (the formula function) is looking for 'z' in the
GlobalEnv, the caller of
> >> apply(), and not finding it. It finds 'y' here explicitly
in the
> >> BoxCoxLambda environment.
> >>
> >> Cheers,
> >> Bert
> >>
> >>
> >>
> >> On Sat, Jul 8, 2023 at 4:28?AM Ron Crump via R-help <r-help at
r-project.org>
> >> wrote:
> >>
> >>> Hi,
> >>>
> >>> Firstly, apologies as I have posted this on
community.rstudio.com too.
> >>>
> >>> I want to optimise a Box-Cox transformation on columns of a
matrix (ie, a
> >>> unique lambda for each column). So I wrote a function that
includes the
> >>> call to MASS::boxcox in order that it can be applied to each
column easily.
> >>> Except that I'm getting an error when calling the
function. If I just
> >>> extract a column of the matrix and run the code not in the
function, it
> >>> works. If I call the function either with an extracted column
(ie dd1 in
> >>> the reprex below) or in a call to apply I get an error (see
the reprex
> >>> below).
> >>>
> >>> I'm sure I'm doing something silly, but I can't
see what it is. Any help
> >>> appreciated.
> >>>
> >>> library(MASS)
> >>>
> >>> # Find optimised Lambda for Boc-Cox transformation
> >>> BoxCoxLambda <- function(z){
> >>> b <- boxcox(lm(z+1 ~ 1), lambda = seq(-5, 5,
length.out = 61), plotit
> >>> = FALSE)
> >>> b$x[which.max(b$y)] # best lambda
> >>> }
> >>>
> >>> mrow <- 500
> >>> mcol <- 2
> >>> set.seed(12345)
> >>> dd <- matrix(rgamma(mrow*mcol, shape = 2, scale = 5), nrow
= mrow, ncol > >>> mcol)
> >>>
> >>> # Try it not using the BoxCoxLambda function:
> >>> dd1 <- dd[,1] # 1st column of dd
> >>> bb <- boxcox(lm(dd1+1 ~ 1), lambda = seq(-5, 5, length.out
= 101), plotit
> >>> = FALSE)
> >>> print(paste0("1st column's lambda is ",
bb$x[which.max(bb$y)]))
> >>> #> [1] "1st column's lambda is 0.2"
> >>>
> >>> # Calculate lambda for each column of dd
> >>> lambdas <- apply(dd, 2, BoxCoxLambda, simplify = TRUE)
> >>> #> Error in eval(predvars, data, env): object 'z'
not found
> >>>
> >>> Created on 2023-07-08 with reprex v2.0.2
> >>>
> >>> Thanks for your time and help.
> >>>
> >>> Ron
> >>> [[alternative HTML version deleted]]
> >>>
> >>> ______________________________________________
> >>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and
more, see
> >>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>> PLEASE do read the posting guide
> >>> http://www.R-project.org/posting-guide.html
> >>> and provide commented, minimal, self-contained, reproducible
code.
> >>>
> >>
> >
> > [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>