Are signifcane asterisks reliable when using robust linear models (rlm) in combination with estimate plots (plot_model)?

[eimichae85]

Hello, due to some outliers I have to work with robust linear models (using the rlm() function, MASS package). I would then like to plot my model results with the plot_model function (package: sjPlot) as I find it a neat way to visualize results.

The summary of the rlm()-based model only provides t-values but not p-values (or significance asterisks). However, when I subject the rlm()-model to the plot_model function, plot_model() somehow provides significance asterisks. 1) How is this possible? 2) Are these significance levels (asterisks indicating p-value ranges) reliable? if so why would I not get them in the summary output of the rlm() model?

Thank you very much for your feedback Michael

Below is an example code and the data to run it

`library(sjPlot) library(MASS)

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

modRob<-rlm(C50~Treatment*Site,data=tab,psi=psi.huber,method="M",maxit=50) summary(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) `

[strengejacke]

See discussion here: https://stat.ethz.ch/pipermail/r-help/2006-July/108659.html

I'd say, as long as you're transparent about the approximation method (and therefore, you could try parameters::model_parameters() and look at the output, which should indicate the approximation method), it's ok to compute them. Compare to mixed models, where you have many different methods (wald, satterthwaite, kenward-rogers, ...), none of which one would say it's perfectly accurate - they're all approximations.

[eimichae85]

Thanks a lot for your swift reply. I was not aware of the possibility to use parameters::model_parameters(model) It works perfectly and gives me indeed p-values. Also the approximation method is mentioned: "p-values (two-tailed) computed using a Wald t-distribution approximation."

I will mention this assuming that this is also what the plot_model() output displays.

Cheers!

[strengejacke]

Yes, plot_model() uses the parameters package to compute the summary tables