2
$\begingroup$

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: enter image description here

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: enter image description here

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)
)

enter image description here

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
$\endgroup$
9
  • $\begingroup$ I'm not sure you provide enough information for any useful suggestions. All we know is that you have 164 observations but want to fit a relatively complex model with 4 fixed and 2 random effects. See this for example: bbolker.github.io/mixedmodels-misc/glmmFAQ.html#gamma-glmms. Starting with a simpler model first (such as Gaussian for log Y) is a reasonable approach. $\endgroup$
    – dipetkov
    Commented Mar 25, 2023 at 22:34
  • $\begingroup$ Thank you for the reply - that link was useful. Eventhough I had looked at the FAQ before, I had missed this detailed section. - I now managed to fit a gamma regression with log link. The residuals, however, look very much like in the plot above with the gaussian and log(response). Is this still acceptable for a model fit? The fact that the data is skewed this way stems from a biological process, and it's just the way it is. $\endgroup$
    – BRB
    Commented Mar 26, 2023 at 13:41
  • $\begingroup$ I also updated my question above now specifying more details. $\endgroup$
    – BRB
    Commented Mar 26, 2023 at 13:53
  • 1
    $\begingroup$ Doesn't the normal model fit better than the gamma? Another idea (a bit out of the left field?) is to use proportional odds (ordinal) regression. You can read more about it here: Ordinal Logistic Regression. It uses the rank of the Y measurements (instead of the raw measurements). You have the extra complexity of multiple measurements per subjects and I'm not sure rms::orm can handle that. $\endgroup$
    – dipetkov
    Commented Mar 26, 2023 at 14:29
  • 1
    $\begingroup$ @dipetkov I think that orm models can handle repeated measurements on individuals, in a marginal modeling context, by starting with a working independence assumption and generating a robust variance-covariance matrix via the rms::robcov() function. Whether that's applicable here could depend on the details of the "group" and "pair" grouping variables, which isn't clear to me from the question. Frank Harrell provides some links to longitudinal ordinal models here. $\endgroup$
    – EdM
    Commented Mar 26, 2023 at 18:36

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Thank you so much for your thoughful explanations, @EdM. "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?" This is a very interesting though. How would I best do that? By trying to fit, e.g. a quadratic relationship? I did indeed expect a linear relationship, but don't no if it really is. Trial, however, is just a code for repeated measures and I do not expect strong changes there, but wanted to still allow for some variation. $\endgroup$
    – BRB
    Commented Mar 27, 2023 at 13:49
  • 1
    $\begingroup$ @BRB You can treat numInfluencers as a categorical predictor (factor) with 5 levels instead of continuous. That's the most flexible. You can respect its ordered nature by coding it as an ordered factor; see the ordered argument in the R help page on factor() and this UCLA web page on the associated orthogonal contrasts that allow for flexibility in the fit. $\endgroup$
    – EdM
    Commented Mar 27, 2023 at 14:37
  • 1
    $\begingroup$ @BRB you could also keep numInfluencers as continuous (including in its interaction term) but add in 3 extra dummy variables indicating whether the actual value is 1 or 2 or 3. With 5 levels you have 4 degrees of freedom to play with, so this gives you 1 linear and 3 categorical predictors. A combined test on the 3 categorical coefficients would evaluate linearity. Section 2.7.3 of Frank Harrell's Regression Modeling Strategies (very useful in general) points to methods for penalized fits including Bayesian shrinkage. $\endgroup$
    – EdM
    Commented Mar 27, 2023 at 14:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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