Aha. Many thanks, John. Never would have gotten there on my own. -- Bert
On Sat, Jul 8, 2023 at 2:01 PM John Fox <j...@mcmaster.ca> wrote: > > Hi Bert, > > On 2023-07-08 3:42 p.m., Bert Gunter wrote: > > Caution: This email may have originated from outside the organization. > > Please exercise additional caution with any links and attachments. > > > > > > 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 > > As it turns out, it's the update() step in boxcox.lm() that fails, and > the update takes place because $y is missing from the lm object, so the > following works: > > BoxCoxLambda <- function(z){ > b <- boxcox(lm(z + 1 ~ 1, y=TRUE), > lambda = seq(-5, 5, length.out = 101), > plotit = FALSE) > b$x[which.max(b$y)] > } > > > > > The identical lambdas do not seem right to me; > > I think that's just an accident of the example (using the BoxCoxLambda() > above): > > > apply(dd, 2, BoxCoxLambda, simplify = TRUE) > [1] 0.2 0.2 > > > dd[, 2] <- dd[, 2]^3 > > apply(dd, 2, BoxCoxLambda, simplify = TRUE) > [1] 0.2 0.1 > > Best, > John > > > 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 <j...@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.4...@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@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@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@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. > >> ______________________________________________ R-help@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.