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.
