QuasiGLM

Adjust Poisson and Binomial Generalised Linear Models to their quasi equivalents for dispersed data

You can install QuasiGLM.jl from the Julia Registry via:

using Pkg
Pkg.add("QuasiGLM")

R has an excellent interface for specifying generalised linear models (GLM) and its base functionality includes a wide variety of probability distributions and link functions. GLM.jl in Julia is also excellent, and boasts a similar interface to its R counterpart. However, in GLM.jl, two key model types are not readily available:

1. quasipoisson
2. quasibinomial

While neither defines an explicit probability distribution, these models are useful in a variety of contexts as they enable the modelling of overdispersion in data. If the data is indeed overdispersed, the estimated dispersion parameter will be >1. Failure to estimate and adjust for this dispersion may lead to inappropriate statistical inference.

QuasiGLM.jl is a simple package that provides intuitive one-line-of-code adjustments to existing Poisson and Binomial GLM.jl models to convert them to their quasi equivalents. It achieves this through estimating the dispersion parameter and using this to make adjustments to standard errors. These adjustments then flow through to updated test statistics, p-values, and confidence intervals.

Here's a Poisson to quasipoisson conversion using the Dobson (1990) Page 93: Randomized Controlled Trial dataset (as presented in the GLM.jl documentation).

using DataFrames, CategoricalArrays, GLM, QuasiGLM

dobson = DataFrame(Counts = [18,17,15,20,10,20,25,13,12], Outcome = categorical([1,2,3,1,2,3,1,2,3]), Treatment = categorical([1,1,1,2,2,2,3,3,3]))

gm = fit(GeneralizedLinearModel, @formula(Counts ~ Outcome + Treatment), dobson, Poisson())
testOutputs = AdjustQuasiGLM(gm, dobson; level=0.95)

And here's a binomial to quasibinomial example using the leaf blotch dataset (McCullagh and Nelder (1989, Ch. 9.2.4)) as seen in multiple textbooks and the SAS documentation:

using DataFrames, CategoricalArrays, GLM, QuasiGLM
    
blotchData = DataFrame(blotch = [0.05,0.00,1.25,2.50,5.50,1.00,5.00,5.00,17.50,0.00,0.05,1.25,0.50,1.00,5.00,0.10,10.00,25.00,0.00,0.05,2.50,0.01,6.00,5.00,5.00,5.00,42.50,0.10,0.30,16.60,3.00,1.10,5.00,5.00,5.00,50.00,0.25,0.75,2.50,2.50,2.50,5.00,50.00,25.00,37.50,0.05,0.30,2.50,0.01,8.00,5.00,10.00,75.00,95.00,0.50,3.00,0.00,25.00,16.50,10.00,50.00,50.00,62.50,1.30,7.50,20.00,55.00,29.50,5.00,25.00,75.00,95.00,1.50,1.00,37.50,5.00,20.00,50.00,50.00,75.00,95.00,1.50,12.70,26.25,40.00,43.50,75.00,75.00,75.00,95.00], variety = categorical(repeat([1,2,3,4,5,6,7,8,9], inner=1, outer=10)), site = categorical(repeat([1,2,3,4,5,6,7,8,9,10], inner=9, outer=1)))
    
blotchData.blotch = blotchData.blotch ./ 100
gm2 = fit(GeneralizedLinearModel, @formula(blotch ~ variety + site), blotchData, Binomial())
testOutputs2 = AdjustQuasiGLM(gm2, blotchData; level=0.95)

Note that results do not exactly equal the R equivalent of GLMs fit with quasibinomial or quasipoisson families. While explorations are continuing, the discrepancy is believed to be the result of differences in optimisation methods in the GLM machinery and floating point calculations.

For example, in the quasipoisson example presented above, the dispersion parameter returned by QuasiGLM.jl and R's glm function with quasipoisson family are equivalent, and the numerical values for the Intercept and Outcome in the summary coefficient table are also equivalent. However, the Treatment variable exhibits different coefficient estimates despite exhibiting the same standard error and p-values.

Here is the R code to test it:

dobson <- data.frame(Counts = c(18,17,15,20,10,20,25,13,12), Outcome = as.factor(c(1,2,3,1,2,3,1,2,3)), Treatment = as.factor(c(1,1,1,2,2,2,3,3,3)))

mod <- glm(Counts ~ Outcome + Treatment, dobson, family = quasipoisson)
summary(mod)

If you use QuasiGLM.jl in your work, please cite it using the following (included as BibTeX file in the package folder):

@Manual{QuasiGLM.jl,
  title={{QuasiGLM.jl}},
  author={Henderson, Trent},
  year={2022},
  month={2},
  url={https://zenodo.org/badge/latestdoi/459915317},
  doi={10.5281/zenodo.6164674}
}