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.