mfasiolo / qgam

Additive quantile regression R package
http://mfasiolo.github.io/qgam/
30 stars 7 forks source link

Stability for low n / low unique predictor values #37

Closed florianhartig closed 4 years ago

florianhartig commented 4 years ago

Hi Matteo,

in DHARMa, I still have problems with the stability of qgam, for low n, or low unique values of the predictors, e.g. https://github.com/florianhartig/DHARMa/issues/172

I'm not sure if the problem is in the quantiles, the gam, or the combination. The only thing I can say is that the package I used before (quantreg) was more stable in these cases (although it may be that the results were not very sensible).

In any case, maybe you can have an eye towards that during development, and / or I wonder if anything could be done. Do you think it would help / make sense to increase the penalty when the model crashes?

Best, Florian

mfasiolo commented 4 years ago

Hi Florian,

Sorry for the delay. If I run this code:

devtools::install_github(repo = "florianhartig/EcoData", subdir = "EcoData", 
                         dependencies = T, build_vignettes = T)
library(EcoData)
library(DHARMa)

fit = glm(feeding ~ attractiveness, birdfeeding, family = "poisson")
summary(fit)
res = simulateResiduals(fit)

I do not get any error. Can you send me an example of what goes wrong with qgam?

Best,

Matteo

florianhartig commented 4 years ago

Sorry, you also have to add

plot(res)

to get the error.

mfasiolo commented 4 years ago

Ok, the problem might come from here:

qgam(res ~ s(pred), data = datTemp, qu = quantiles[i])

which gives the error

Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) : 
  A term has fewer unique covariate combinations than specified maximum degrees of freedom

This come from mgcv::gam, and signal the fact that datTemp$pred takes only 5 unique values (in this case), while s(pred) uses 10 knots by default. The easiest thing to do would be to modify your code by doing for instance:

qgam(res ~ s(pred, k = min(length(unique(datTemp$pred)), 10)), data = datTemp, qu = quantiles[i])
florianhartig commented 4 years ago

Hi Matteo,

excellent, this makes sense, I will change this!

Cheers, Florian

florianhartig commented 4 years ago

Works perfect, many thanks!