Help fitting glmm with positive, right skewed, continuous data

Question

I am trying to fit a glmm in R, with a right-skewed response variable that is theoretically continuous, but in my case ranging between 0.4 and 1.8 with more lower values (it's a biological measurement).

I also want to include 3 categorical and one integer as predictors, and have two random grouping variables as well. This is the data structure:

group               pair           trial       var            treatment           type              numInfluencers 
 Length:164         Length:164         1:42   Min.   :0.4472   Length:164         Length:164         Min.   :0.000  
 Class :character   Class :character   2:56   1st Qu.:0.4876   Class :character   Class :character   1st Qu.:1.000  
 Mode  :character   Mode  :character   3:66   Median :0.5538   Mode  :character   Mode  :character   Median :2.000  
                                              Mean   :0.6630                                         Mean   :2.085  
                                              3rd Qu.:0.7049                                         3rd Qu.:3.000  
                                              Max.   :1.7882                                         Max.   :4.000  

and a histogram of my response:

I tried fitting it with gamma distribution, but just keep getting lots of warnings (In (function (start, objective, gradient = NULL, hessian = NULL, ... : NA/NaN function evaluation), presumably meaning that the model is not converging, I guess.

The only somewhat reasonable fit so far I achieved with a gaussian model and log transforming the response (using a log-link instead was considerably worse). But transforming the variable would not be my preferred approach, and the residual are also still looking less than ideal:

Any advice on how to approach this and what else I could try? Thank you.

Update:

I now fit a gamma model with log link with a residual plot looking very much like the one above.

m <- glmmTMB(var  ~  treatment * numInfluencers + type + trial + 
               (1|pair)+(1|group),
             data = df,family=Gamma(link=log))

Interestingly, I get a very different residual plot when I fit the same model with lme4:

m <- glmer(var  ~  treatment * numInfluencers + type + trial + 
               (1|pair)+(1|group),
             data = df,family=Gamma(link=log)
)

Why is that? And is this type of residual plot even meaningful here? A QQ plot of this model does not look so bad, actually.

Lastly, I am aware that ideally there would be more data for this model. If I try running this model anyway, what would be the best way of assessing the goodness of fit and decide whether the model is "reasonable"?

This is the model output that I am getting right now. I am slightly suspicious because all the variables are significant, but this is also possible. I tried adding a non meaningful variable into it and it did not show up as significant.

Family: Gamma  ( log )
Formula:          var ~ treatment * numInfluencers + type + trial + (1 | pair) +      (1 | group)
Data: df

     AIC      BIC   logLik deviance df.resid 
   -94.5    -63.5     57.3   -114.5      154 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 pair   (Intercept) 0.01799  0.1341  
 group  (Intercept) 0.02580  0.1606  
Number of obs: 164, groups:  pair, 46; group, 4

Dispersion estimate for Gamma family (sigma^2): 0.0565 

Conditional model:
                                Estimate Std. Error z value Pr(>|z|)   
(Intercept)                     -0.10100    0.10008  -1.009  0.31286   
treatmenttreated                -0.23177    0.09579  -2.420  0.01554 * 
numInfluencers                  -0.05665    0.01931  -2.934  0.00335 **
typeB                           -0.15956    0.05516  -2.893  0.00382 **
trial2                          -0.17589    0.05727  -3.071  0.00213 **
trial3                          -0.18381    0.05972  -3.078  0.00209 **
treatmenttreated:numInfluencers  0.08041    0.03922   2.050  0.04035 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Answer 1

I suspect that the problems with your residual plots come from some mis-specification of the model. The DHARMa vignette shows some ways to look for such problems via its simulated residuals.

In particular, the numInfluencers has integer values that range from 0 to 4. You are treating that as having an exactly linear association with outcome at treatment=0 (the individual coefficient for numInfluencers) and a linear interaction with treatment. You might need to fit that more flexibly: is the effect of moving from 0 to 1 influencers really expected to be the same as moving from 3 to 4? Also, the trial number in your model has a fixed association with outcome independent of treatment; I wonder whether that's realistic.

With respect to the choice between the log-transformed and Gamma models, as a biologist I personally find the log transformations easier to think about. That might be because I've been thinking about log transformations for about 60 years, ever since I learned the relationship between pH and hydrogen-ion activity, but I've not had much cause to think about how to interpret a Gamma-distributed variable with a log link between its mean and a linear predictor. If your audience understands that, or you can explain it simply to them, then the Gamma model does get around the potential confusion with the log-transformed data, where you are modeling the mean of the logs instead of the log of the mean.

My comment to @dipetkov about ordinal regression was more for her benefit, as you seem to be handling the repeated measures well enough with your mixed models. The orm() function doesn't handle random effects. Another way to deal with repeated measures, without modeling random effects explicitly, is to fit a model as if the repeated measures don't matter (a "working independence" assumption), then calculate a "robust" or "sandwich" variance-covariance matrix for the coefficient estimates based on the "clusters" in your data, that takes the repeated measures into account.

You might apply that in your situation if you treated group as a fixed effect and treated the pairs as the clustering variable. That's what the robcov() function can do with an orm model. (With only 4 levels of group you might consider treating it as fixed in your mixed model, anyway.) I think that the Bayesian ordinal mixed modeling in the rmsb package only allows for a single random effect, so you might have to treat group as a fixed effect were you to try that.