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
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 therms::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$