Great, the answer linked is a great explanation. Thank you! So basically if I conduct an explorative analysis and I do not know whether Z interacts with X or it is its mediator I would need to run independent regressions to understand that. In both cases, if I do not include the Z variables (either as interaction or as mediator) in an initial model, I will be in both cases able to capture a smaller fraction of Y variance. So adding them either as Y:X or as M will help me to "model" Y better. Do I understand this correctly?

Thank you very much. I think what I do not quite get is whether these concepts are mutually exclusive. In other words, I can have a variable X that influences Y, influences (partially) M and then Y, and also affects Y by its interaction with M. Am I seeing this completely wrong or there is a way to tease this apart? In mediation, do we assume that M is completely caused by X? Thank you

Absolutely, true, I was just trying to write the link functions myself, and I realized this makes no sense at all! Thank you very much for looking into this.

I guess my problem is that I do not quite get how to implement the back transformation similarly to what Ben Bolker did in the link above

Thanks, I just saw it, quite a complex task. I would need o basically write a link function myself and then parse into tran. I guess something like: stackoverflow.com/questions/15931403/…

thank you, the point is that I need to account for the geometric mean of the data. I might create my own function but then I am not sure how to provide this to tran.

Yes, you are right, what I wrote does not match the data! Sorry about that

Thank you. I think I kind of got the solution. First I define the interested factors e.g. A=c(0,0,0,0,1,0) and B=c(0,1,0,0,0,0,0). Then I call contrast(emm1, method = list("A-B" = A - B)). This seemed to work