strengejacke / sjPlot

sjPlot - Data Visualization for Statistics in Social Science
https://strengejacke.github.io/sjPlot
606 stars 90 forks source link

Backtransformation of sqrt() transformed estimates using plot_model() #943

Open eimichae85 opened 2 weeks ago

eimichae85 commented 2 weeks ago

Hello, I would like to use the plot_model() command to plot the estimates & SE of the following model

modRob<-lm(sqrt(C50)~Treatment*Site,data=tab)

I would like to get the untransformed/backtransformed estimates (not the sqrt-transformed estimates) in my final plot_model output. How do I do that? Under help(plot_model) it says: "Estimates of linear models remain untransformed. Use transform= NULL if you want the raw, non-transformed estimates". However, if I use "transform=NULL" I still get the sqrt-transformed values. Also "transform="^2" does not work.

What command do I have to use to back-transform my sqrt-transformed estimates?

Thanks a lot for your help Michael

Below is the code and a dummy data set I am working with

#DATA
tab<-structure(list(Site = structure(c(2L, 1L, 1L, 1L, 2L, 2L, 1L, 
                                       2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 
                                       2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 
                                       1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 
                                       1L), levels = c("Alpthal", "Birmensdorf"), class = "factor", contrasts = "contr.sum"), 
                    Treatment = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                            2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                            3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
                                            4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                            1L, 1L), levels = c("control", "4000_100", "4000_30", "4000_50"
                                            ), class = "factor"), C50 = c(490139.1569, 210684.3876, 177220.7609, 
                                                                          144855.0927, 143589.8017, 346457.7683, 208220.2951, 763146.691, 
                                                                          337140.0098, 473545.559, 133118.6341, 196323.6835, 324372.0377, 
                                                                          660669.6405, 221500.1217, 130437.8918, 228093.4435, 64947.9752, 
                                                                          98456.18456, 108942.3415, 258502.6415, 38607.72435, 60942.19915, 
                                                                          276165.4695, 193098.7689, 215195.3215, 175787.9858, 200215.9341, 
                                                                          80513.00079, 120383.9266, 127098.0298, 0, 36507.1922, 63337.42105, 
                                                                          119012.168, 136744.0301, 243138.5785, 237136.1503, 143608.6254, 
                                                                          196767.0193, 72837.33976, 37430.74448, 180981.6339, 21085.30347, 
                                                                          62652.19874, 430516.4452, 101975.2863, 286981.5731, 174188.796, 
                                                                          148825.2492, 376042.6975, 0, 361023.8245, 73044.75764, 317343.18, 
                                                                          155962.495)), row.names = c(NA, -56L), class = "data.frame")

##Define the factors
tab$Site<-as.factor(tab$Site)
tab$Treatment<-as.factor(tab$Treatment)

contrasts(tab$Site) <- "contr.sum"
contrasts(tab$Treatment)<- "contr.sum"

tab$Treatment <- relevel(tab$Treatment, ref="control") #only use for graphics, NOT statistics

library(sjPlot)
library(MASS)

modRob<-lm(sqrt(C50)~Treatment*Site,data=tab)
summary(modRob)

#parameters::model_parameters(modRob)

plot_model(modRob, vline.color = "red", colors = c("black","black"),
           sort.est=F,show.values = TRUE,
           value.offset = .45,value.size = 6,title = "",
           dot.size=5.3, line.size =0.8 ,
           transform = NULL , order.terms = c(2,3,1,4,6,7,5), auto.label = FALSE)
eimichae85 commented 2 weeks ago

Hi, Even after a further search I could not figure this out. Is anybody able to help me solve my issue above (i.e. what command should I use for "transform=...." to backgransform the estimates after a sqrt transformation "^2" does not work)?

Thanks a lot Michael

strengejacke commented 2 weeks ago

Can you try

power_of_2 <- function(x) x^2

plot_model(modRob,
  vline.color = "red", colors = c("black", "black"),
  sort.est = F, show.values = TRUE,
  value.offset = .45, value.size = 6, title = "",
  dot.size = 5.3, line.size = 0.8,
  transform = "power_of_2",
  order.terms = c(2, 3, 1, 4, 6, 7, 5), auto.label = FALSE
)
eimichae85 commented 2 weeks ago

Hi thanks a ton for your reply.

I already tried this approach. Unfortunately, if I use your code the estmiate-values and the corresponding SEs are not aligned (I think the power-transformation is only applied to the SE bars). See image below

Capture