diff --git a/R/description-glm.R b/R/description-glm.R
index 3c157ea..91eb0cf 100644
--- a/R/description-glm.R
+++ b/R/description-glm.R
@@ -6,6 +6,7 @@ descript.GLMModelResult <- function(x, ...){
   } else {
     etype <- "g-computation-type"
   }
+
   if(is.character(x$settings$g_family)){
     if(x$settings$g_family == "nb"){
       family <- "negative binomial with unknown dispersion"
@@ -20,7 +21,16 @@ descript.GLMModelResult <- function(x, ...){
 
   if(!is.null(x$data$formula)){
 
-    form <- x$mod$formula
+    if(is.character(x$settings$g_family)){
+      if(x$settings$g_family == "nb"){
+        form <- x$data$formula
+      } else {
+        form <- x$mod$formula
+      }
+    } else {
+      form <- x$mod$formula
+    }
+
     output <- c(
       output,
       sprintf("Treatment group mean estimates from a GLM working model of family %s and link %s using formula: \n",
diff --git a/tests/testthat/test-glm.R b/tests/testthat/test-glm.R
index af67183..19f14cc 100644
--- a/tests/testthat/test-glm.R
+++ b/tests/testthat/test-glm.R
@@ -241,6 +241,21 @@ test_that("GLM full function -- NEGATIVE binomial, permuted block", {
     covariate_to_include_strata=TRUE)
   expect_equal(class(non), "GLMModelResult")
 
+  non
+
+  # Test with formula
+  non <- robincar_glm2(
+    df=DATA2,
+    response_col="y",
+    treat_col="A",
+    car_strata_cols=c("z1"),
+    formula="y ~ A + z1",
+    car_scheme="permuted-block",
+    g_family="nb",
+    g_accuracy=7)
+
+  non
+
 })
 
 test_that("GLM Settings", {