Combining 'make.tran' transformations?
a-hurst opened this issue · 4 comments
Hi!
First of all, sincere thanks for creating and maintaining this wonderful package. It makes my workflow massively easier and it's one of the first packages I tell people about who are new to stats in R.
For my research I do a lot of analyses involving reaction time data, which means that I typically log-transform responses before analysis. Additionally, because I run models with brms and it makes it easier to set priors, I also z-score the log-transformed RTs prior to running the model. emmeans makes it easy to work with these models in their transformed log + scale units, but because I've done two different transformations on the data before modelling it's not as simple as using make.tran to get estimates back in milliseconds.
I've figured out how to get it to work with a rough hack (overwriting the linkinv function in the list returned by make.tran with a custom function that re-scales and exponentiates the values), but I'm wondering if there's a simpler/cleaner way to do multiple transformations that I've missed? If not, would you consider adding an extra argument to make.tran("scale", ...) that allows specifying an additional simple transformation (e.g. "log", "sqrt") that was performed prior to scaling and adds a corresponding exp() or ** 2 for the inverse function?
Thanks in advance!
- Austin
It turns out that this is fairly easy to do using recursive coding tactics...
> y = rnorm(100, mean = 15)
> tran1 = make.tran("genlog")
> tran = make.tran("scale", y = tran1$linkfun(y), inner = tran1)
> tran
$linkfun
\(mu) outer$linkfun(inner$linkfun(mu))
<bytecode: 0x0000019989b413f8>
<environment: 0x000001998b8fe9f8>
$linkinv
\(eta) inner$linkinv(outer$linkinv(eta))
<bytecode: 0x0000019989b35b08>
<environment: 0x000001998b8fe9f8>
$mu.eta
\(eta) inner$mu.eta(outer$linkinv(eta))*outer$mu.eta(eta)
<bytecode: 0x0000019989b462f8>
<environment: 0x000001998b8fe9f8>
$valideta
\(eta) inner$valideta(outer$linkinv(eta))
<bytecode: 0x0000019989b382c0>
<environment: 0x000001998b8fe9f8>
$name
[1] "scale(2.77, 0.053)(log(mu + 1))"
I need to write documentation for inner, and also test it out -- especially to see whether "scale" works right, as it requires some special handling. Then I'll push it up.
BTW, a side effect is that make.tran() now also accepts arguments for stats::make.link() e.g. make.tran("sqrt") and a bunch of others e.g. make.tran("atanh")
OK, I think I have it all working now
> tran = make.tran("scale", y = warpbreaks$breaks, inner = "log")
> mod = with(tran,
+ lm(linkfun(breaks) ~ wool*tension, data = warpbreaks))
> emmeans(mod, ~ tension|wool)
wool = A:
tension emmean SE df lower.CL upper.CL
L 1.0909 0.285 48 0.517 1.665
M -0.2852 0.285 48 -0.859 0.289
H -0.2832 0.285 48 -0.857 0.291
wool = B:
tension emmean SE df lower.CL upper.CL
L 0.0939 0.285 48 -0.480 0.668
M 0.1556 0.285 48 -0.418 0.729
H -0.7719 0.285 48 -1.346 -0.198
Results are given on the scale(3.24, 0.437)(log) (not the response) scale.
Confidence level used: 0.95
> emmeans(mod, ~ tension|wool, type = "response")
wool = A:
tension response SE df lower.CL upper.CL
L 41.2 5.13 48 32.0 52.9
M 22.6 2.81 48 17.6 29.0
H 22.6 2.82 48 17.6 29.0
wool = B:
tension response SE df lower.CL upper.CL
L 26.6 3.32 48 20.7 34.2
M 27.4 3.41 48 21.3 35.2
H 18.2 2.28 48 14.2 23.4
Confidence level used: 0.95
Intervals are back-transformed from the scale(3.24, 0.437)(log) scale
Please note that I changed this a bit from my previous example, in how y is handled. When we are using "scale" with "inner", I have it set to replace y with inner$linkfun(y) in the recursive call, so that we are talking about the right compound transformation of y.
@rvlenth Looks great, thanks for taking a look at this so quickly! Will definitely make my workflow easier.
Closing this issue as it is resolved, and the new version that supports this is now on CRAN