Search code examples
rtransformationnlmeemmeans

Specifying that model is logit transformed to plot backtransformed trends


I have fitted a lme model in R with a logit transformed response. I have not been able to find a direct command that does the logit transformation so I have done it manually.

logitr<-log(r/1-r)

I then use this as response in my lme model with interaction between two factors and a numerical variable.

model<-lme(logitr<-factor1*factor2*numeric,random=1|random)

Now, R obviously do not know that this model is logit transformed. How can I specify this to R?

I have without luck tried:

update(model, tran="logit")

The reason why I want to specify that the model is logit transformed is because I want to plot the backtransformed results using the function emmip in the emmeans package, showing the trends of the interaction between my variables.

Normally (if I only had factors) I would just use:

update_refgrid_model<-update(ref_grid(model, tran="logit"))

But this approach does not work when I want to use emmip to plot the trends of the interaction between a numerical variable and factors. If I specify:

emmip(update_refgrid_model, factor1~numeric|factor2, cov.reduce = range, type = "response")

then I do not get any trends plotted, only the estimate for the average level on the numerical variable.

So, how can I specify the logit transformation and plot the backtransformed trends of a lme model with factors interacting with numerical variables?


Solution

  • You don't update the model object, you update the reference grid:

    rg = update(ref_grid(model, cov.reduce = range), tran = "logit")
    emmip(rg, factor1~numeric|factor2, type = "response")
    

    It is possible to update a model with other things, just not the transformation; that is in the update method for emmGrid objects.

    Update

    Here's an example showing how it works

    require(emmeans)
    ## Loading required package: emmeans
    
    foo = transform(fiber, p = (strength - 25)/25)
    foo.lm = lm(log(p/(1-p)) ~ machine*diameter, data = foo)
    
    emm = emmeans(foo.lm, ~diameter|machine, 
                  tran = "logit", at = list(diameter = 15:32))
    ## Warning in ref_grid(object, ...): There are unevaluated constants in the response formula
    ## Auto-detection of the response transformation may be incorrect
    
    emmip(emm, machine ~ diameter)
    

    emmip(emm, machine ~ diameter, type = "r")
    

    Created on 2020-06-02 by the reprex package (v0.3.0)