# How to interpret significant factors in a GLMM gamma model that appears to be good but predicts poorly the observations?

I have a data set of reaction times (resp.lat) of 120 subjects (su) and 36 observations per subject, which I am fitting with a GLMM model with gamma distribution and identity link as recommended by Lo & Andrew (2015) (the identity link can create problems with gamma but it does not appear to be the case here). The model is

resp.lat ~ stake.i.c + dir.c + win.prob.c + (stake.i.c + dir.c + win.prob.c || su)

The predictors are continuous variables taking a small number of possible values (two for dir.c, five for win.prob.c and five for stake.i.c). I am potentially interested in identifying random effects to correlate them with other subject data at a later stage of analysis.

Update: The predictors were originally categorical variables but exploratory plots showed that a linear relation was a reasonable. While this coding scheme assumes linear relationships between RTs and predictors, it uses less parameter than categorical variables. On exploratory plots (above), the effect of the variables is dwarfed by the response variability, which makes the significance of the fixed effects all the more surprising (below). I have removed previous information referring to the categorical variables.

Update 2 I have added DHARMa Residual plots as suggested by Dimitris Rizopoulos (see below).

Looking at the summary, the model appears to fit well the data:

> fitRespLat11 <- glmer(resp.lat ~ stake.i.c + dir.c + win.prob.c +
(stake.i.c + dir.c + win.prob.c || su),
data=tmp, family = Gamma(link = "identity"),
control=glmerControl(optimizer =  "bobyqa"))
> summary(fitRespLat11)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma  ( identity )
Formula: resp.lat ~ stake.i.c + dir.c + win.prob.c + (stake.i.c + dir.c + win.prob.c || su)
Data: tmp
Control: glmerControl(optimizer = "bobyqa")

AIC      BIC   logLik deviance df.resid
58026.3  58083.4 -29004.2  58008.3     4160

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.7918 -0.6890 -0.2162  0.4502  4.8900

Random effects:
Groups   Name        Variance  Std.Dev.
su       (Intercept) 5.413e+03  73.576
su.1     stake.i.c   2.008e+03  44.806
su.2     dir.c       2.944e+03  54.260
su.3     win.prob.c  7.514e+04 274.108
Residual             2.172e-01   0.466
Number of obs: 4169, groups:  su, 120

Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept)  612.589      6.290  97.390  < 2e-16 ***
stake.i.c    -61.037      6.001 -10.171  < 2e-16 ***
dir.c        -33.530      5.481  -6.117 9.52e-10 ***
win.prob.c    34.071      8.262   4.124 3.72e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) stk..c dir.c
stake.i.c  -0.089
dir.c      -0.014  0.126
win.prob.c  0.191 -0.226 -0.108


The model fit without problems.

> isSingular(fitRespLat11)
[1] FALSE
> rePCA(fitRespLat11)
$su Standard deviations (1, .., p=4): [1] 588.19493 157.88241 116.43287 96.14736 Rotation (n x k) = (4 x 4): [,1] [,2] [,3] [,4] [1,] 0 1 0 0 [2,] 0 0 0 1 [3,] 0 0 1 0 [4,] 1 0 0 0 attr(,"class") [1] "prcomplist"  The model fits without warnings and random effects appear appropriate: > anova(fitRespLat11,fitRespLat12,fitRespLat13,fitRespLat14) Data: tmp Models: fitRespLat14: resp.lat ~ stake.i.c + dir.c + win.prob.c + (1 | su) fitRespLat13: resp.lat ~ stake.i.c + dir.c + win.prob.c + (stake.i.c || su) fitRespLat12: resp.lat ~ stake.i.c + dir.c + win.prob.c + (stake.i.c + dir.c || su) fitRespLat11: resp.lat ~ stake.i.c + dir.c + win.prob.c + (stake.i.c + dir.c + win.prob.c || su) Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq) fitRespLat14 6 58203 58241 -29096 58191 fitRespLat13 7 58129 58173 -29057 58115 76.699 1 < 2.2e-16 *** fitRespLat12 8 58056 58107 -29020 58040 74.820 1 < 2.2e-16 *** fitRespLat11 9 58026 58083 -29004 58008 31.496 1 1.999e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  LRT and BIC suggest that all random effects are statistically significant. In summary, the model converged without warnings, it is not singular and all random effects are statistically significant (LRTs). The residual errors appears to be small (SD=0.466) and all fixed effects are statistically significant (t values > 4, Wald test < 0.0001). However, these "results" go against what the plot show. The plot of the predicted values against the observations shows a lot of unexplained variability. plot(fitRespLat11,resp.lat ~ fitted(.))  The DHARMa diagnostic plot (residuals were simulated with use.u=TRUE to avoid negative values with gamma identity link) also shows a bad fit, which I interpret as showing strong underdispersion simOut <- simulateResiduals(fittedModel = fit.glmer , n = 250, use.u=TRUE) plot(simOut)  I have several (slightly updated) questions: 1. Why are the results statistically significant ? Dimitris Rizopoulos mentioned in the comment that it might be possible when the fit is bad but I don't understand what is happening intuitively in this case. 2. How should I change the model so that I get expected results ? Why is gamma not working. When I look at the estimated dispersion, using a gamma disribution does not seem a bad choice a priori. 3. I am surprised by the fact that the residual variance in the model is very low (SD=0.466; see below), which could explain why factors are statistically significant. If the SD is on a different scale than the random effects as it seems, how do I transform it? I would expect something like sd(tmp$resp.lat - fitted(fitRespLat11))
286.0733

Histogram of reaction times fitted GLM with and gamma distribution

fit.glm <- glm(resp.lat ~  1 , data=dat, family = Gamma(link = "identity"))
mu <- coef(fit.glm)
disp <- summary(fit.glm)$$dispersion hist(dat$$resp.lat,n=100,xlim=c(0,2000),freq=FALSE,main="",xlab="RT ms")
lines(density(dat$$resp.lat)) x<-seq(0,max(dat$$resp.lat),length=100)
y <- dgamma(x, 1/disp, scale = mu * disp)
lines(x,y,col=2)


Underdispersion suggest that I have less variability than I should. Indeed, the density plot (black line) is less wide than the fitted gamma distribution (red line).

Response time distribution distribution for each subject (ordered according to mean RT):

The gamma dispersion parameters estimated by glm and glmer are similar

 > summary(fit.glm)\$dispersion
[1] 0.2628206
> disp  <- sigma(fitRespLat11)^2
[1] 0.2171703


• Thank you. I will have a look on DHARMa. I have tried to find functions to compute R2 (MuMIn::r.squaredGLMM and piecewiseSEM::rsquared, rsquared) but they cannot compute the variance for GLMM gamma with identity link. The simple function r2.corr.mer in GLMM Faq yields 22%, which seems reasonable given the above figure). – gavril Feb 6 '19 at 15:32