Robust options to fit GLM or GAM for overdispersed Poisson counts (quasipoisson or negative binomial) in R

Question

It appears glmrob from the robustbase library does not support quasipoisson or 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 or 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.

Answer 1

It's not too hard to make your own quasi-adjustment function, which can be layered on top of a robust estimation procedure (or any procedure that returns Wald estimates for the standard errors of the parameters and from which we can extract residuals that are suitable for computing an estimate of dispersion).

• extract the coefficient table
• compute the dispersion, $\phi$, usually as the ratio of the sum of squared Pearson residuals to the residual df (number of obs - number of coefficients)
• adjust the standard errors by multiplying by $\sqrt{\phi}$
• recalculate the $Z$ (or $t$) scores and the $p$-values on the basis of the adjusted standard errors

The code below hould be fairly general: the main warning is that not all model types allow you to extract coefficient tables via coef(summary(.)), and the column names might not be consistent across model types (using broom::tidy might be safer)

quasi_tab <- function(fit) {
    df.resid <- nobs(fit) - length(coef(fit))
    dispersion <- sum(residuals(fit, type = "pearson")^2)/df.resid
    cc <- coef(summary(fit))
    cc[,"Std. Error"] <-     cc[,"Std. Error"]*sqrt(dispersion)
    cc[,"z value"] <- cc[,"Estimate"]/cc[,"Std. Error"]
    ## calculate two-sided p-values (2x lower tail)
    cc[,"Pr(>|z|)"] <- 2*pnorm(-1*abs(cc[,"z value"]), lower.tail = TRUE)
    return(cc)
}

printCoefmat(quasi_tab(fit))

Answer 2

It is important to distinguish between

1. standard errors that are robust to mis-specifying the variance (which is usually "fixed" by argument robust=T

2. a fitted model where all distribution parameters, (e.g. mean, variance) are robust to outliers in the data set. [Aeberhard et al. (2021) look at a robust gamlss fit, (although I don't know if it includes the negative binomial, and I have not tested it myself so I don't know how well it works.]