[R-sig-ME] [RE] Response time modeling: Warning message with non-normal RAs

Ulf Köther ukoether @ending from uke@de
Wed Oct 17 18:28:30 CEST 2018


Dear Gabriel,

as no one of the real professionals stepped in yet, I will try to
provide some thoughts:

(1) In my experience, Gamma GLMMs with an identity link are very buggy
all the time, so normally, when I am confronted with positive continuous
data, I mostly never consider them in the first place...why not use
another link function when using GLMMs, because that's exactly the goal
they were invented for: modeling data that does not belong to a normal
distribution by using an appropriate link function.

(2) I see why you don't want to transform the raw RT data and use an
ordinary LMM, but in my opinion, using a Gamma GLMM with a log-Link is a
valid option. Especially because you can insert the raw RT as a DV and
then model it on the scale of the linear predictor and after that, get
all parameter estimates and (which is much more valuable in my opinion)
the predicted values on the scale of the dependent variable without any
hard work, just use the invert link function (e.g., use the effects
package with lme4 and get the raw RT predictions). The interpretation if
such a model is much easier as you can directly get the predictions on
the scale of the DV.

Is it totally reasonable with regards to power? That's another question
that a real RT specialist should answer...

(2) I know the article you mentioned and skimmed it again. I really
wonder why they do not mention a log-link function as a valid option and
just compare the LMM with the transformed DV against a Gamma GLMM with
the identity link or an inverse Gaussian with identity link. It seems to
my that the Gamma GLMM with a log link is never really discussed. But,
as they also use the inverse Gaussian with the following link function
in lme4:

invfn <- function() {
    ## link
    linkfun <- function(y) -1000/y
    ## inverse link
    linkinv <- function(eta)  -1000/eta
    ## derivative of invlink wrt eta
    mu.eta <- function(eta) { 1000/(eta^2) }
    valideta <- function(eta) TRUE
    link <- "-1000/y"
    structure(list(linkfun = linkfun, linkinv = linkinv,
                   mu.eta = mu.eta, valideta = valideta,
                   name = link),
              class = "link-glm")
}

Why not try this when you want to go along the line the authors suggest?
The parameter estimates in their analyses seem really comparable between
Gamma and inverse Gaussian GLMMs to me. Just have a look at the
supplement! ;-)

Try to play around with your generated data with the log-link function
and with the inverse Gaussian to see which one recovers your parameters
best and go from there...!

Greetings,

Ulf




Am 17.10.2018 um 11:06 schrieb Baud-Bovy Gabriel:
> My apologies for sending this post again.  I am not sure that it went
> through because I did not
> receive it myself (though I think my options should have let me receive
> it).  Also, there was
> a line missing due to hurried cut-and-paste in the example below.
> 
> Dear all,
> 
> I have a dataset with response times of 200 participants measured 36
> times in different conditions.
> The RTs distribution is right skewed as it is generally the case for
> RTs. I would like to avoid transform them
> (see Lo & Andrew
> (2015)](https://www.frontiersin.org/articles/10.3389/fpsyg.2015.01171/full).
> However, I often encounter warning messages with this dataset when using
> a GLMM with Gamma
> distribution and identity link. The problem is not caused by an over
> parametrization as it arises
> even in the case in a intercept only model.
> 
>  > fit <- glmer(choice.lat ~ 1 + (1|su), data=dat,
> family=Gamma(link="identity"))
> Warning messages:
> 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>    Model failed to converge with max|grad| = 0.220808 (tol = 0.001,
> component 1)
> 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>    Model is nearly unidentifiable: very large eigenvalue
>   - Rescale variables?
> 
>  > fit
> Generalized linear mixed model fit by maximum likelihood (Laplace
>    Approximation) [glmerMod]
>   Family: Gamma  ( identity )
> Formula: choice.lat ~ 1 + (1 | su)
>     Data: dat
>        AIC       BIC    logLik  deviance  df.resid
>   22629.26  22649.78 -11311.63  22623.26      6914
> Random effects:
>   Groups   Name        Std.Dev.
>   su       (Intercept) 0.8477
>   Residual             0.5721
> Number of obs: 6917, groups:  su, 200
> Fixed Effects:
> (Intercept)
>        3.044
> convergence code 0; 2 optimizer warnings; 0 lme4 warnings
> 
> The dispersion parameter of the estimated Gamma distribution is
> sigma(fit.glmer)^2 = 0.327
> with a population mean (Fixed effect) of circa  3.0 sec which fits
> reasonably well the RTs.
> 
> After some poking around, I figured out that the problem probably comes
> from the
> fact that the participants means are NOT distributed normally.  In fact,
> their distribution
> is similarly skewed. Removing the two extreme subjects or using a log
> link function "solves" the convergence
> problem/ warning messages but I don't like either solution because the
> is no valid reason to exclude these subjects
> and using a log transform cause theoretical problems when one want to
> interpret
> effects of more complex models (see Lo & Andrew (2015)).
> 
> I have seen many discussions of about warning message/convergence problem (
> e.g.
> https://rstudio-pubs-static.s3.amazonaws.com/33653_57fc7b8e5d484c909b615d8633c01d51.html).
> The usual advice is to try with different optimizers and/or simplify the
> model, which do not help here (they all give
> similar warning messages; see below).  I have two questions:
> 
> 1) Can somebody explain intuitively why non-normal random effects lead
> to these warning
>     messages ?  I would have naively assumed that it would simply lead
> to large a large variance estimate.
> 
> 2) Is there something within the lme4 framework  to tackle this problem?
> Should I simply ignore the
>      warning messages ?
> 
> 3) What principle approach would you suggest ? I have  seen that brms
> allow in principle to specify
>     non-normal random effect (I think that only student  is implemented
> so far) but I haven't tried yet.
> 
> Here is a toy problem that reproduces the issue (it depends on the seed
> actually). Depending
> on the choice of the dispersion parameter for the REs, it is possible
> also to induce convergence
> problem.
> 
> # artificial data set
> nsu <- 200
> nrep <- 36
> mu <- 3.0      # mean RT
> phi.su <- 0.8  # dispersion participants
> phi <- 0.3     # dispersion
> d <- expand.grid(
>    su = paste0("S",1:nsu),
>    rep = 1:nrep,
>    mu =NA,
>    rt = NA)
> 
> set.seed(125)
> # generate 200 Gamma distributed  participant means
> mu.su <- sort(rGamma(nsu, shape=1/phi.su, scale=phi.su*mu))
> d$mu <- mu.su[d$su]
> # add Gamma noise
> d$rt <- rGamma(nsu*nrep,shape=1/phi,scale=phi*d$mu)
> 
> # fit intercept GLMM model (gives warning messages)
> fit <- glmer(rt ~  1 + (1|su), data=d, family=Gamma(link="identity"))
> 
> # histogram of participant means and response times
> par(mfrow=c(1,2))
> hist(mu.su,n=50)
> hist(d$rt,n=200,xlim=c(0,10))
> 
>  > all_fit(fit)
> bobyqa. : [OK]
> Nelder_Mead. : [OK]
> optimx.nlminb : [OK]
> optimx.L-BFGS-B : [OK]
> nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
> nloptwrap.NLOPT_LN_BOBYQA : [OK]
> nmkbw. : [OK]
> ....
> were 13 warnings (use warnings() to see them)
> 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>    Model failed to converge with max|grad| = 0.215844 (tol = 0.001,
> component 1)
> 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>    Model is nearly unidentifiable: very large eigenvalue
>   - Rescale variables?
> 3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>    Model failed to converge with max|grad| = 0.220808 (tol = 0.001,
> component 1)
> 4: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>    Model is nearly unidentifiable: very large eigenvalue
>   - Rescale variables?
> 5: In optwrap(optimizer, devfun, start, rho$lower, control = control,  ... :
>    convergence code 1 from optimx
> 6: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>    Model failed to converge with max|grad| = 0.215713 (tol = 0.001,
> component 1)
> 7: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>    Model is nearly unidentifiable: very large eigenvalue
>   - Rescale variables?
> 8: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>    Model failed to converge with max|grad| = 0.102562 (tol = 0.001,
> component 1)
> 9: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>    Model failed to converge with max|grad| = 0.109779 (tol = 0.001,
> component 1)
> 10: In checkConv(attr(opt, "derivs"), opt$par, ctrl =
> control$checkConv,  ... :
>    Model failed to converge with max|grad| = 0.109779 (tol = 0.001,
> component 1)
> 11: In optwrap(optimizer, devfun, start, rho$lower, control = control,
> ... :
>    convergence code 2 from nmkbw
> 12: In checkConv(attr(opt, "derivs"), opt$par, ctrl =
> control$checkConv,  ... :
>    Model failed to converge with max|grad| = 0.220811 (tol = 0.001,
> component 1)
> 13: In checkConv(attr(opt, "derivs"), opt$par, ctrl =
> control$checkConv,  ... :
>    Model is nearly unidentifiable: very large eigenvalue
>   - Rescale variables?
> 
> 
> Thank you,
> 
> Gabriel
> 
> --
> ---------------------------------------------------------------------
> Gabriel Baud-Bovy               tel.: (+39) 348 172 4045     (mobile)
> UHSR University                       (+39) 02 2643 3429 (laboratory)
> via Olgettina, 58                     (+39) 02 9175 1540  (secretary)
> 20132 Milan, Italy              email:       gabriel.baud-bovy using hsr.it
> ---------------------------------------------------------------------
> 
> 
> 
> Rispetta l’ambiente: non stampare questa mail se non è necessario.
> Respect the environment: print this email only if necessary.
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
--

_____________________________________________________________________

Universitätsklinikum Hamburg-Eppendorf; Körperschaft des öffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard Göke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Prölß, Martina Saurin (komm.)
_____________________________________________________________________

SAVE PAPER - THINK BEFORE PRINTING


More information about the R-sig-mixed-models mailing list