glmrob from the
robustbase library does not support
negative binomial for the moment for the robust fitting of overdispersed Poisson counts. The same is true for
robust::glmRob. Is anyone aware of good options in R to fit GLMs or GAMs to overdispersed Poisson counts using a robust method? Maybe
gamlss could be an option, https://link.springer.com/article/10.1007/s11222-020-09979-x#Sec16, via the
robust=TRUE option in
gamlss, https://rdrr.io/cran/GJRM/man/gamlss.html? Or any other options, ideally with an example? I also need the capability to add an offset term to the model formula (this is also not currently supported by
glmRob) and I need to be able to calculate confidence and prediction intervals on model predictions for custom covariate values (for prediction intervals resampling coefficients from a multivariate normal distribution based on the coefficients & variance covariance matrix and then using the percentile method, ie using population prediction intervals, would be OK for me). Just using the
sandwich package to obtain Huber-White standard errors is also not quite enough for me, as that would just give better estimates of the standard errors, but not influence my actual parameter estimates (when in fact what I want is that my parameters would not be influenced too much by outliers in the data).
The specific problem I am working on is the modelling of mortality data, using a periodic
splines2::mSpline term to model seasonality. Occasionally there is years that have more severe flu waves in winter, causing higher mortality then & resulting in biases in the coefficients if one fits on a limited nr of years. This gets better if I use a robust Poisson GLM. But then I was worried that I wasn't taking into account overdispersion.
A reproducible example below:
# weekly mortality by sex & age from https://osf.io/k84rz/ library(archive) load(archive_read("https://osf.io/download/ydpvu/", file = 3)) # MOCY database # dataframe mocy mocy$age_group = factor(mocy$age_start, labels=unique(attr(mocy$age_start, "names"))) mocy$sex = factor(mocy$sex) head(mocy) country = "Belgium" region = "Belgium" data = mocy[mocy$country_name==country&mocy$region_name==region,]
A negative binomial fit to this mortality data would take into account overdispersion but is not robust to outliers in the data :
library(MASS) library(splines2) fit_negbin_from2010 = glm.nb(deaths ~ mSpline(week, df=5, Boundary.knots=c(0,52.17857), periodic=TRUE) + # period M-spline to model seasonality (age_group+sex+year)^2 + offset(log(personweeks)), data=data[data$year>=2010&data$year<2020,])
A robust Poisson GLM would be insensitive to outliers in the data (and hence could be fit on all years, even including the pandemic years, where we had large mortality outliers), but does not explicitly model overdispersion:
library(robustbase) fit_robpois_from2010_allyears = glmrob(deaths ~ mSpline(week, df=5, Boundary.knots=c(0,52.17857), periodic=TRUE) + # periodic M-spline to model seasonality (age_group+sex+year)^2 + log(personweeks), # offset term not supported by glmrob, so putting it in as a term instead family=poisson(log), method="Mqle", weights.on.x="hat", data=data[data$year>=2010,], control=glmrobMqle.control(maxit=1000))
So any good approach that would be robust plus explicitly take into account overdispersion??
EDIT: Stefan Van Aelst just sent me this article (forthcoming in the Journal of the American Statistical Association): https://www.dropbox.com/s/r831vu18ruy4r4y/RDE_accepted.pdf?dl=1, which seems to be close to what I need. Together with this R package: https://github.com/RdeGlmLassoGam/RdeGlmLassoGam/blob/main/man/RDE.Rd and some examples here https://github.com/RdeGlmLassoGam/RDE_data_examples/blob/main/Data_Analyse/Example_ili.visits.R... And a negative binomial GAMLSS https://rdrr.io/cran/GJRM/man/gamlss.html may also be a good solution with option
robust=TRUE, using method as in https://link.springer.com/article/10.1007/s11222-020-09979-x.
glmnb()? (BTW, questions only about R are often considered off topic.) $\endgroup$