9
$\begingroup$

I have long-term collection data, and I'd like to test, whether the number of animals collected is influenced by weather effects. My model looks like below:

glmer(SumOfCatch ~ I(pc.act.1^2) +I(pc.act.2^2) + I(pc.may.1^2) + I(pc.may.2^2) + 
                   SampSize + as.factor(samp.prog) + (1|year/month), 
      control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=1e9,npt=5)), 
      family="poisson", data=a2)

Explanation of the used variables:

  • SumOfCatch: number of animals collected
  • pc.act.1, pc.act.2: axes of a principal component representing weather conditions during sampling
  • pc.may.1, pc.may.2: axes of a PC representing weather conditions in May
  • SampSize: number of pitfall traps, or collecting transects of standard lengths
  • samp.prog: method of sampling
  • year: year of sampling (from 1993 to 2002)
  • month: month of sampling (from Aug to Nov)

The fitted model's residuals show considerable inhomogeneity (heteroscedasticity?) when plotted against fitted values (see Fig.1):

Residuals vs. Fitted values - Main model

My main question is: is this a problem making the reliability of my model questionable? If so, what can I do to resolve it?

So far I have tried the followings:

  • control for overdispersion by defining observation-level random effects, i.e. using a unique ID for each observation, and applying this ID variable as random effect; although my data do show considerable overdispersion, this did not help as the residuals became even more ugly (see Fig. 2)

Residuals vs. Fitted values - Model with OD control

  • I fitted models without random effects, with quasi-Poisson glm and glm.nb; also yielded similar residual vs. fitted plots to the original model

As far as I know, there might be ways for the estimation of heteroscedasticity-consistent standard errors, but I have failed to find any such method for Poisson (or any other kind of) GLMMs in R.

In response to @FlorianHartig: the number of observations in my dataset is N=554, I think this is a fair obs. number for such a model, but of course, the more the merrier. I post two figures, first of which is the DHARMa scaled residual plot (suggested by Florian) of the main model.

enter image description here

The second figure is from a second model, in which the only difference is that it contains the observation-level random effect (the first does not).

enter image description here

UPDATE

Figure of the relationship between a weather-variable (as predictor, i.e. x-axis) and sampling success (response):

Weather-PC and sampling success

UPDATE II.

Figures showing predictor values vs. residuals:

Predictors vs. Residuals

Improve this question
$\endgroup$
7
  • $\begingroup$ Have you considered running a nonparametric estimator? Or comparing an ols to a median regression? I realize that poisson is the dominant model in bio but GLMs are inconsistent under heteroscedasticity and OLS is not. $\endgroup$
    – Superpronker
    Dec 15, 2016 at 7:31
  • 1
    $\begingroup$ Sometimes overdispersion is caused by zero-inflation. In that case you could try a poisson model with zero-inflation parameter or a hurdle model. The glmmADMB package has great features to deal with this: glmmadmb.r-forge.r-project.org/glmmADMB.html $\endgroup$
    – Niek
    Dec 15, 2016 at 9:00
  • $\begingroup$ Dear @Superpronker thanks for the suggestion, I did not checked OLS, I did not realize that this approach would be flexible enough to handle my data. I'll look into it $\endgroup$
    – zozi9126
    Dec 15, 2016 at 11:58
  • $\begingroup$ Dear @Niek in my data, there are no observations of zero - otherwise I did think of zeroinfl and hurdle models (in package 'pscl') due to their good handling of overdispersion, but they are only usable on data with zeros in the response. Back a few months ago I did try glmmADMB, but it did not yield better results. Cheers, ZR $\endgroup$
    – zozi9126
    Dec 15, 2016 at 12:01
  • 1
    $\begingroup$ @mdewey the rationale behind this, is that the relationship between weather effects and sampling success follows an optimum: the probability and extent of sampling success is highest in a given range (in this case, zero and its vicinity) of the predictor. When the predictor's value is further away from this optimum, the sampling success will be lower, corresponging to a suboptimum. I use the quadratic term, because (1) predictors are rescaled and recentered at zero, and (2) this gives a better approximation to a linear connection. Cheers, ZR $\endgroup$
    – zozi9126
    Dec 15, 2016 at 13:51

1 Answer 1

Reset to default
10
$\begingroup$

It is difficult to assess the fit of the Poisson (or any other integer-valued GLM for that matter) with Pearson or deviance residuals, because also a perfectly fitting Poisson GLMM will exhibit inhomogeneous deviance residuals.

This is especially so if you do GLMMs with observation-level REs, because the dispersion created by OL-REs is not considered by the Pearson residuals.

To demonstrate the issue, the following code creates overdispersed Poisson data, that is then fitted with a perfect model. The Pearson residuals look very much like your plot - hence, it may be that there is no problem at all.

This problem is solved by the DHARMa R package, which simulates from the fitted model to transform the residuals of any GL(M)M into a standardized space. Once this is done, you can visually assess / test residual problems, such as deviations from the distribution, residual dependency on a predictor, heteroskedasticity or autocorrelation in the normal way. See the package vignette for worked-through examples. You can see in the lower plot that the same model now looks fine, as it should.

If you still see heteroscedasticity after plotting with DHARMa, you will have to model dispersion as a function of something, which is not a big problem, but would likely require you to move to JAGs or another Bayesian software. 

library(DHARMa)
library(lme4)

testData = createData(sampleSize = 200, overdispersion = 1, randomEffectVariance = 1, family = poisson())

fittedModel <- glmer(observedResponse ~ Environment1 + (1|group) + (1|ID), family = "poisson", data = testData, control=glmerControl(optCtrl=list(maxfun=20000) ))

# standard Pearson residuals
plot(fittedModel, resid(., type = "pearson") ~ fitted(.) , abline = 0)

# DHARMa residuals
plot(simulateResiduals(fittedModel))

enter image description here

enter image description here

Improve this answer
$\endgroup$
12
  • $\begingroup$ Dear @FlorianHartig! Thank You for Your insights, I tried plotting with DHARMa. Based on the plot there is still something up, causeing the lower quantile to be shaped like a reciprocal curve, rather than a straight line. You have mentioned that in this case, a solution might be to model dispersion as a function of something - could You help in exactly how can i assess such a function? Cheers, ZR $\endgroup$
    – zozi9126
    Dec 15, 2016 at 11:52
  • $\begingroup$ Can you send me or post the plot? Some minor variability is expected when your N is small $\endgroup$
    – Florian Hartig
    Dec 15, 2016 at 13:12
  • $\begingroup$ Dear @FlorianHartig the question was edited, now showing the DHARMa plots as well! $\endgroup$
    – zozi9126
    Dec 15, 2016 at 15:39
  • $\begingroup$ @Z.Radai - what I see in the plots is that your residuals are systematically too high for low model predictions. This looks more like a problem of model structure (missing predictors?) than a distribution problem - I would try plotting the residuals against possible and potentially missing predictors. $\endgroup$
    – Florian Hartig
    Dec 15, 2016 at 16:43
  • 2
    $\begingroup$ I would not be concerned about heteroskedasticity, in your case it's moderate and the effect on the inference should be mild - the only issue I see is the systematic underprediction for small values, which will not be solved by modeling the variance. But if you must know see here stats.stackexchange.com/questions/247183/… $\endgroup$
    – Florian Hartig
    Dec 21, 2016 at 9:26

Your Answer

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

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