hendersontrent · hendersontrent · Feb 17, 2022 · Feb 17, 2022
diff --git a/README.md b/README.md
@@ -13,29 +13,55 @@ Pkg.add("QuasiGLM")
 
 ## Motivation
 
-`R` has an excellent interface for specifying [generalised linear models](https://en.wikipedia.org/wiki/Generalized_linear_model) (GLM) and its base functionality includes a wide variety of probability distributions and link functions. [`GLM.jl`](https://juliastats.org/GLM.jl/v0.11/) in `Julia` is also excellent, and boasts a similar (and also excellent!) interface. However, in `GLM.jl`, two key models are not readily available:
+`R` has an excellent interface for specifying [generalised linear models](https://en.wikipedia.org/wiki/Generalized_linear_model) (GLM) and its base functionality includes a wide variety of probability distributions and link functions. [`GLM.jl`](https://juliastats.org/GLM.jl/v0.11/) 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 inaccurate statistical inference.
+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. It achieves this through adjustments to standard errors which flows through to updated test statistics, *p*-values, and confidence intervals.
+`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.
 
 ## Usage
 
-Adjustments can be made to existing `GLM` models in `Julia` simply using `QuasiGLM.jl`:
+Here's a Poisson to quasipoisson conversion using the Dobson (1990) Page 93: Randomized Controlled Trial dataset (as presented in the [`GLM.jl` documentation](https://juliastats.org/GLM.jl/v0.11/#Fitting-GLM-models-1)).
 
 ```
 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](https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glimmix_sect016.htm):
+
+```
+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)
+```
+
+### Comparison to R results
+
+Note that results do not exactly equal the `R` equivalent of GLMs fit with `quasibinomial` or `quasipoisson` families. While explorations are continuing, it is believed to be the result of differences in the 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, `Treatment` variable exhibits different `Estimate` values 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)
+```
+
 ## Citation instructions
 
 If you use `QuasiGLM.jl` in your work, please cite it using the following (included as BibTeX file in the package folder):

diff --git a/src/AdjustQuasiGLM.jl b/src/AdjustQuasiGLM.jl
@@ -1,8 +1,6 @@
 """
     AdjustQuasiGLM(model, ϕ; level)
-
 Estimates dispersion parameter, adjusts original GLM to reflect the dispersion and returns results in a pretty DataFrame.
-
 Usage:
 ```julia-repl
 AdjustQuasiGLM(model, ϕ; level)

diff --git a/src/PeripheralFunctions.jl b/src/PeripheralFunctions.jl
@@ -1,8 +1,6 @@
 """
     PearsonResiduals(model, data)
-
 Calculates Pearson residuals between model predicted values and the actual response variable values.
-
 Usage:
 ```julia-repl
 PearsonResiduals(model, data)
@@ -32,9 +30,7 @@ end
 
 """
     EstimateDispersionParameter(residuals, model)
-
 Estimates the dispersion parameter ϕ by standardising the sum of squared Pearson residuals against the residual degrees of freedom.
-
 Usage:
 ```julia-repl
 EstimateDispersionParameter(residuals, model)
@@ -56,9 +52,7 @@ end
 
 """
     coefarray(model, ϕ; level)
-
 Calculates relevant statistics for inference based of estimates and dispersion-adjusted standard errors and returns results in a concatenated array.
-
 Usage:
 ```julia-repl
 coefarray(model, ϕ; level)
@@ -80,4 +74,4 @@ function coefarray(model::StatsModels.TableRegressionModel, ϕ::Real; level::Rea
     ci = se * quantile(TDist(dof_residual(model)), (1 - level) / 2)
     ct = hcat(coefnames(model), cc, se, tt, p, cc + ci, cc - ci)
     return ct
-end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -20,4 +20,23 @@ using DataFrames, CategoricalArrays, GLM, Distributions, PrettyTables, QuasiGLM,
 
     testOutputs = AdjustQuasiGLM(gm, dobson; level=0.95)
     @test testOutputs isa DataFrames.DataFrame
+
+    #--------------
+    # Quasibinomial
+    #--------------
+
+    # Set up data and divide percentage by 100 to get proportion
+
+    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
+
+    # Fit binomial model
+
+    gm2 = fit(GeneralizedLinearModel, @formula(blotch ~ variety + site), blotchData, Binomial())
+
+    # Correct standard errors using quasi correction
+
+    testOutputs2 = AdjustQuasiGLM(gm2, blotchData; level=0.95)
+    @test testOutputs2 isa DataFrames.DataFrame
 end