3
$\begingroup$

I've built a gam using mgcv and the following code

m1 <- gam(y ~ s(x1, k = 10) + s(x2, k = 55), data = df, method = "REML",
            family = poisson(link = "sqrt"),
            select = TRUE)

y is discrete species counts at a given location and is right skewed, and x1 and x2 represent some environmental correlates.

The model summary seems fine:

Family: poisson 
Link function: sqrt 

Formula:
y ~ s(x1, k = 10) + s(x2, k = 55)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 15.16921    0.01001    1516   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
             edf Ref.df Chi.sq p-value    
s(x1)  8.989      9  62237  <2e-16 ***
s(x2) 52.185     54  32089  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.531   Deviance explained = 57.9%
-REML =  78296  Scale est. = 1         n = 2497

And the diagnostic plots also seem fine, with the exception of the qqplot as it doesn't follow the red line, and the scale on the y-axis is an order of magnitude greater than that on the x-axis

enter image description here

How can I interpret the qqplot with regards to the model interpretation? Have I specified an incorrect family or link function, or is the model not performing as well as summary(m1) and the other plots would have me believe?

|cite|improve this question
$\endgroup$
2
$\begingroup$

The QQ-plot suggests overdispersion, but the residuals vs linear predictor plot suggests that you don't have the correct model specification, or at least you should check this.

Is there a reason for using the square-root link function? Is it better with the log-link, which is the more standard link function in this setting?

If the problem persists after switching to the log-link, I'd see if things improve by assuming the response is conditionally distributed Negative Binomial, and try fitting the model with family = nblink = 'log').

Beyond that I think you'll need to start plotting residuals against your covariates and consider if you have missing covariates or structure that you didn't tell the model about (are the samples clustered or grouped in some way?). I'd also check that the basis dimensions for your smooths were large enough; the EDF of the smooths in the output from summary() is quite close to the maximum possible given the values of k you have stated - there's not been much penalization here.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ That is an expected residuals vs linear predictor plot for a Poisson GLM: the variance is equal to the mean. $\endgroup$ – AdamO Nov 8 '18 at 17:12
  • $\begingroup$ specifying family = nb(link = "log") solved the issues, the qqplot (more or less) now looks how it should and the residuals v predictor plot looks like a random scatter. I have checked the model with gam.check() and all the edf values for the smooths are less than k' with non-significant p-values $\endgroup$ – KaanKaant Nov 8 '18 at 21:50
  • $\begingroup$ @AdamO my understanding is that, for what is shown in that top-right plot (deviance residuals plotted against $\eta_i$), we might expect roughly constant spread if we have the right conditional distribution for the response. That we don't see that suggests that we have the wrong mean-variance relationship. I didn't mean to say that this is was sufficient to conclude that we have the right/wrong model; at low count data sets we may see non-constant variance in this plot even if we do have the correct model. $\endgroup$ – Gavin Simpson Nov 8 '18 at 22:23
1
$\begingroup$

A QQ plot makes no sense for a Poisson GLM using standard normal distribution for theoretical quantiles. The theoretical quantiles for Poisson errors are related to the variance stabilizing transform $\sqrt{Y}$. See here for how to inspect the distributional assumptions of a GLM with deviance residuals.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ But that QQ-plot doesn't use the standard normal distribution to generate theoretical quantiles (well, not with the options the user selected - if they'd used method = 'normal' then the standard normal would have been used for the theoretical quantiles). What is being done here is the direct randomisation method of Augustin et al (2012). The data simulation method from the same paper is also available via method = 'simulate'. There's more on this in ?mgcv::qq.gam too. $\endgroup$ – Gavin Simpson Nov 8 '18 at 22:04
  • $\begingroup$ (And I only know this because I wrote the function(s) [not shown by the OP] that were used to produce the plots they showed.) Perhaps the "Theoretical quantiles" title for the x-axis is misleading in this case, but in my defence I stuck with what Simon Wood labelled the axis in qq.gam() to avoid the usual user questions of why this was different - it didn't occur to me until now that this might be a source of confusion too! $\endgroup$ – Gavin Simpson Nov 8 '18 at 22:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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