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 <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.