rvlenth / emmeans

Estimated marginal means
https://rvlenth.github.io/emmeans/
358 stars 31 forks source link

Combining 'make.tran' transformations? #462

Closed a-hurst closed 6 months ago

a-hurst commented 8 months ago

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

rvlenth commented 8 months ago

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")

rvlenth commented 8 months ago

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.

a-hurst commented 8 months ago

@rvlenth Looks great, thanks for taking a look at this so quickly! Will definitely make my workflow easier.

rvlenth commented 6 months ago

Closing this issue as it is resolved, and the new version that supports this is now on CRAN