6
$\begingroup$

I would like some clarification whether my model is well specified or not (since I do not have much experience with Beta regression models).

My variable is the percentual of dirth area in the denture. For every pacient, the dentist applied a special product in either left or right side on the denture (leaving the other side as placebo) in order to remove dirth area.

After that, he calculate the total area of each side of the denture, and the total dirth area for each side.

I need to test whether the product is efficient to remove the dirth.

My initial model (prop.bio is the proportion of dirth area):

library(glmmTMB)    
m1 <- glmmTMB(prop.bio ~ Product*Side + (1|Pacients), data, family=list(family="beta",link="logit"))

Update:

My final model after manual backward selection via TRV test (and it is also the main question of the researcher):

m1.f <- glmmTMB(prop.bio ~ Product + (1|Pacients), data, family=list(family="beta",link="logit"))

My residual diagnosis using DHARMa:

library(DHARMa)
res = simulateResiduals(m1.f)
plot(res, rank = T)

enter image description here

According to my reading on DHARMa vignette, my model could be wrong based on the right plot. What should I do then? (Is my model specification wrong?)

Thanks in advance!

Data: 

structure(list(Pacients = structure(c(5L, 6L, 2L, 11L, 26L, 29L, 
20L, 24L, 8L, 14L, 19L, 7L, 13L, 4L, 3L, 5L, 6L, 2L, 11L, 26L, 
29L, 20L, 24L, 8L, 14L, 19L, 7L, 13L, 4L, 3L, 23L, 25L, 12L, 
21L, 10L, 22L, 18L, 27L, 15L, 9L, 17L, 28L, 1L, 16L, 23L, 25L, 
12L, 21L, 10L, 22L, 18L, 27L, 15L, 9L, 17L, 28L, 1L, 16L), .Label = c("Adlf", 
"Alda", "ClrW", "ClsB", "CrCl", "ElnL", "Gema", "Héli", "Inác", 
"Inlv", "InsS", "Ircm", "Ivnr", "Lnld", "Lrds", "LusB", "Mart", 
"Mrnz", "Murl", "NGc1", "NGc2", "Nlcd", "Norc", "Oliv", "Ramr", 
"Slng", "Svrs", "Vldm", "Vlsn"), class = "factor"), Area = c(3942, 
3912, 4270, 4583, 2406, 2652, 2371, 4885, 3704, 3500, 4269, 3743, 
3414, 4231, 3089, 4214, 3612, 4459, 4678, 2810, 2490, 2577, 4264, 
4287, 3487, 4547, 3663, 3199, 3836, 3237, 3846, 4116, 3514, 3616, 
3609, 4053, 3810, 4532, 4380, 4103, 4552, 3745, 3590, 3386, 3998, 
4449, 3367, 3698, 3840, 4457, 3906, 4384, 4000, 4156, 3594, 3258, 
4094, 2796), Side = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L), .Label = c("Right", "Left"), class = "factor"), Biofilme = c(1747, 
1770, 328, 716, 1447, 540, 759, 1328, 2320, 1718, 1226, 977, 
1193, 2038, 1685, 2018, 1682, 416, 679, 2076, 947, 1423, 1661, 
1618, 1916, 1601, 1833, 1050, 1780, 1643, 1130, 2010, 2152, 812, 
2550, 1058, 826, 1526, 2905, 1299, 2289, 1262, 1965, 3016, 1630, 
1823, 1889, 1319, 2678, 1205, 472, 1694, 2161, 1444, 1062, 819, 
2531, 2310), Product = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L), .Label = c("No", "Yes"), class = "factor"), prop.bio = c(0.443176052765094, 
0.452453987730061, 0.0768149882903981, 0.156229543966834, 0.601413133832086, 
0.203619909502262, 0.320118093631379, 0.271852610030706, 0.626349892008639, 
0.490857142857143, 0.287186694776294, 0.261020571733903, 0.349443468072642, 
0.481682817300874, 0.545483975396568, 0.478879924062648, 0.465669988925803, 
0.0932944606413994, 0.145147498931167, 0.738790035587189, 0.380321285140562, 
0.552192471866511, 0.389540337711069, 0.377420107301143, 0.549469457986808, 
0.352100285902793, 0.5004095004095, 0.328227571115974, 0.464025026068822, 
0.507568736484399, 0.293811752470099, 0.488338192419825, 0.612407512805919, 
0.224557522123894, 0.706566916043225, 0.261041204046385, 0.216797900262467, 
0.336716681376876, 0.66324200913242, 0.316597611503778, 0.502855887521968, 
0.3369826435247, 0.547353760445682, 0.890726520968695, 0.407703851925963, 
0.409755001123848, 0.561033561033561, 0.356679286100595, 0.697395833333333, 
0.270361229526587, 0.12083973374296, 0.386405109489051, 0.54025, 
0.347449470644851, 0.295492487479132, 0.251381215469613, 0.618221787982413, 
0.82618025751073)), row.names = c(NA, -58L), class = "data.frame")
Improve this question
$\endgroup$

3 Answers 3

Reset to default
10
$\begingroup$

tl;dr it's reasonable for you to worry, but having looked at a variety of different graphical diagnostics I don't think everything looks pretty much OK. My answer will illustrate a bunch of other ways to look at a glmmTMB fit - more involved/less convenient than DHARMa, but it's good to look at the fit as many different ways as one can.

First let's look at the raw data (which I've called dd):

library(ggplot2); theme_set(theme_bw())
ggplot(dd,aes(Product,prop.bio,colour=Side))+
    geom_line(colour="gray",aes(group=Pacients))+
    geom_point(aes(shape=Side))+
    scale_colour_brewer(palette="Dark2")

enter image description here

My first point is that the right-hand plot made by DHARMa (and in general, all predicted-vs-residual plots) is looking for bias in the model, i.e. patterns where the residuals have systematic patterns with respect to the mean. This should never happen for a model with only categorical predictors (provided it contains all possible interactions of the predictors), because the model has one parameter for every possible fitted value ... we'll see below that it doesn't happen if we look at fitted vs residuals at the population level rather than the individual level ...

The quickest way to get fitted vs residual plots (e.g. analogous to base-R's plot.lm() method or lme4's plot.merMod()) is via broom.mixed::augment() + ggplot:

library(broom.mixed)
aa <- augment(m1.f, data=dd)
gg2 <- (ggplot(aa, aes(.fitted,.resid))
    + geom_line(aes(group=Pacients),colour="gray")
    + geom_point(aes(colour=Side,shape=Product))
    + geom_smooth()
)

enter image description here

These fitted and residual values are at the individual-patient level. They do show a mild trend (which I admittedly don't understand at the moment), but the overall trend doesn't seem large relative to the scatter in the data.

To check that this phenomenon is indeed caused by predictions at the patient rather than the population level, and to test the argument above that population-level effects should have exactly zero trend in the fitted vs. residual plot, we can hack the glmmTMB predictions to construct population-level predictions and residuals (the next release of glmmTMB should make this easier):

aa$.fitted0 <- predict(m1.f, newdata=transform(dd,Pacients=NA),type="response")
aa$.resid0 <- dd$prop.bio-aa$.fitted0
gg3 <- (ggplot(aa, aes(.fitted0,.resid0))
    + geom_line(aes(group=Pacients),colour="gray")
    + geom_point(aes(colour=Side,shape=Product))
    + geom_smooth()
)

(note that if you run this code, you'll get lots of warnings from geom_smooth(), which is unhappy about being run when the predictor variable [i.e., the fitted value] only has two unique levels)

enter image description here

Now the mean value of the residuals is (almost?) exactly zero for both levels (Product=="No" and Product=="Yes").

As long as we're at it, let's check the diagnostics for the random effects:

lme4:::dotplot.ranef.mer(ranef(m1.f)$cond)

enter image description here

This looks OK: no sign of discontinuous jumps (indicating possible multi-modality in random effects) or outlier patients.

other comments

  • I disapprove on general principles of reducing the model based on which terms seem to be important (e.g. dropping Side from the model after running anova()): in general, data-driven model reduction messes up inference.
Improve this answer
$\endgroup$
3
  • $\begingroup$ "in general, data-driven model reduction messes up inference." +1! $\endgroup$
    – Daniel
    Aug 23, 2019 at 6:36
  • $\begingroup$ Tks Ben! I just edited my question saying that I did backward selection via Likelihood ratio test (And in some other problems, I use LASSO/RIDGE). What do you consider being the best approach when I want to simplify my model? $\endgroup$
    – Guilherme Parreira
    Aug 26, 2019 at 18:55
  • 1
    $\begingroup$ If you simplify your model based on the response variable, in any way, you will make inference much more difficult. If you're only interested in prediction, and not in interval estimates/confidence intervals/p-values, that doesn't matter. see e.g. stats.stackexchange.com/questions/377527/… (although there may be better answers floating around on CV) $\endgroup$
    – Ben Bolker
    Aug 27, 2019 at 11:41
5
$\begingroup$

Have a look at the section about glmmTMB in the vignette of DHARMa. It seems to be an issue with regard to how predictions are calculated given the random effects.

As an alternative, you may try the GLMMadaptive package. You can find examples using the DHARMa here.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Yes. I think you should quote this from the message that appears when you run simulateResiduals(): "glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot" $\endgroup$
    – Ben Bolker
    Aug 22, 2019 at 19:51
5
$\begingroup$

I am the developer of DHARMa. Dimitris and Ben are correct, the pattern originates from the known issue that glmmTMB does not (yet) allow making predictions based on fixed effects only, which sometimes produces this pattern. I hope we can fix this issue with the next release of glmmTMB, which should allow fixed-effect predictions.

[EDIT Nov 21: this problem was fixed in glmmTMB approximately 1yr ago.]

In your case, it is obvious that the predicted variable in your model is based on fixed and random effects, because your fixed effects have only one categorical predictor, so you should have only 2 values on your x axis. We can produce a plot using only fixed effects as predictors easily by hand:

plotResiduals(data$Product, res$scaledResiduals)

Which results in a plot that looks fine

enter image description here

btw, agree with Ben that I would not do model selection based on significance, this is essentially p-hacking. If you start with Product*Side, report this model, unless you think there is a serious issue with the inference.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Tks Florian! I just edited my question saying that I did backward selection via Likelihood ratio test (And in some other problems, I use LASSO/RIDGE). Is it still wrong? What should I do? I always look forward a solution to simplify the interpretation for the researcher. $\endgroup$
    – Guilherme Parreira
    Aug 26, 2019 at 18:54
  • 1
    $\begingroup$ Just stay with the initial model. This doesn't appear to be a case where MS is necessary. $\endgroup$
    – Florian Hartig
    Sep 28, 2019 at 18:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.