search

florianhartig / DHARMa

Diagnostics for HierArchical Regession Models
http://florianhartig.github.io/DHARMa/
203 stars 22 forks source link

lme4::simulate.merMod produces Inf / NA for inverse.gaussian(link='identity') or Gamma(link='inverse') families #222

Open BeffaraBee opened 3 years ago

BeffaraBee commented 3 years ago

Hi there,

I made multiple trials to use DHARMa on this GLMM: mod1=glmer(RT~DISPSL+(DISPSL|subj),family=Gamma(link="inverse"),CorrectBehavior)

Then I asked for models checks: testDispersion(mod1) simulationOutput <- simulateResiduals(fittedModel = mod1, plot = T)

... and got this error: Simulations from your fitted model produce infinite values. Consider if this is sensible Error in securityAssertion("Simulations from your fitted model produce NaN values. DHARMa cannot calculated residuals for this. This is nearly certainly an error of the regression package you are using", : Message from DHARMa: During the execution of a DHARMa function, some unexpected conditions occurred. Even if you didn't get an error, your results may not be reliable. Please check with the help if you use the functions as intended. If you think that the error is not on your side, I would be grateful if you could report the problem at https://github.com/florianhartig/DHARMa/issues Context: Simulations from your fitted model produce NaN values. DHARMa cannot calculated residuals for this. This is nearly certainly an error of the regression package you are using In addition: Warning message: In rgamma(nsim * length(ftd), shape = shape, rate = shape/ftd) : NAs produced

What is strange is that it works with the same model, simply using a LMM.

Any idea...?

florianhartig commented 3 years ago

As the error message says, there seems to be something going wrong in the glmer simulations, so this is an lme4, not a DHARMa problem. It's not strange that it works using lmer, as the problem lies probably in the inverse link.

For saying more, I would need your data. Also / alternatively, you could try to run simulate(mod1) to see where things go wrong (although, as I said, this is probably a numeric / programming problem on the lme4 side)

As a quick fix, you could try to run simulate residuals with option re.form = ~0 , which switches off the re-simulation of the REs (see https://www.rdocumentation.org/packages/DHARMa/versions/0.3.3.0/topics/simulateResiduals and https://www.rdocumentation.org/packages/lme4/versions/1.1-26/topics/predict.merMod). I suspect that the problem somehow lies in the combination of the RE simulations and the inverse link.

BeffaraBee commented 3 years ago

Hi,

Thanks for your prompt answer.

DATA MOVED TO https://github.com/florianhartig/DHARMa/blob/master/Code/DHARMaIssues/222.txt

Thanks a lot for your help!

psdmp commented 3 months ago

Hi Fabien and BeffaraBee,

I wondered if you could update this issue with a potential solution or further suggestions.

I receive the same error in the following model that includes an interaction: mod2c <- glmer(n1ballrt ~ t + cond + t*cond + (1|id), family = inverse.gaussian(log), data=d.1b.all.rt, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5))) n1ballrt = reaction time on a 1-bck task t = time (pretest, posttest) cond = condition (active, control) id = participant id

I tried the following distributions:

  1. inverse.gaussian(identity)
  2. inverse.gaussian(log)
  3. Gamma(link = "inverse")
  4. gaussian(identity)

1-#3 result in the same error in DHARMa:

"simulationOutput <- simulateResiduals(fittedModel = mod2c, n = 1000) Simulations from your fitted model produce infinite values. Consider if this is sensible Error in securityAssertion("Simulations from your fitted model produce NaN values. DHARMa cannot calculated residuals for this. This is nearly certainly an error of the regression package you are using", : Message from DHARMa: During the execution of a DHARMa function, some unexpected conditions occurred. Even if you didn't get an error, your results may not be reliable. Please check with the help if you use the functions as intended. If you think that the error is not on your side, I would be grateful if you could report the problem at https://github.com/florianhartig/DHARMa/issues

Context: Simulations from your fitted model produce NaN values. DHARMa cannot calculated residuals for this. This is nearly certainly an error of the regression package you are using."

However, specifying a model without the interaction does not yield this error: mod2b <- glmer(n1ballrt ~ t + cond + (1|id), family = inverse.gaussian(log), data=d.1b.all.rt, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))

Assessing models without the interaction term shows underdispersion when using DHARMa, which is only resolved when using Gaussian distribution with the identity link. This distribution is not appropriate to my model. I enclosed two pictures of the DHARMa plots: 1. for a model without an interaction with an inverse.gaussian distribution with log link function; 2. a model with a Gaussian distribution with identity link function.

I am using package lme4 version 1.1.35.3, R version 4.4.0 and R Studio version 2024.04.1+748

I am a novice in using generalized mixed effects models, and DHARMa package, and hope to learn more while I am using the packages for my data. Any suggestions would be much appreciated and warmly welcome.

Thank you!

1

Screen Shot 2024-06-01 at 10 28 29 PM

2

Screen Shot 2024-06-01 at 10 33 04 PM
florianhartig commented 3 months ago

Hi,

first of a all, apologies to @BeffaraBee - I dropped the ball here, which is that I forgot to follow up on the data you provided.

I have looked that the example of @BeffaraBee - here, the issue was, as I suspected, that if you simulate with the REs unconditional as is the default in DHARMa, you will get into areas where the Gamma distribution is not finite. Calculating conditional simulations solves the issue, although the residuals didn't look good. 

dat = read.table("222.txt", header = T, sep = "|")
head(dat)
library(lme4)
mod1=glmer(RT~DISP+(DISP|subj),family=Gamma(link="inverse"),data = dat)

simulationOutput <- simulateResiduals(fittedModel = mod1, plot = T)
# produces infinity error 

x = predict(mod1,  type = "response")
min(x)
max(x)

# unconditional simulations, default in DHARMa
simulate(mod1)
simulate(mod1, re.form = ~0)

# conditional simulations 
simulate(mod1, re.form = NULL)

simulationOutput <- simulateResiduals(fittedModel = mod1, plot = T, re.form = NULL)

@psdmp, I can't say much about your data unless you provide it, ideally via dput, but it on a very superficial glance, it seems to me that in both cases inappropriate or problematic combinations of link / distribution could be an issue.