[R] Getting an error calling MASS::boxcox in a function
John Fox
j|ox @end|ng |rom mcm@@ter@c@
Sat Jul 8 23:01:27 CEST 2023
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 <jfox using 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 using 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 using 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 using 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 using 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.
>>
More information about the R-help
mailing list