diff --git a/docs/_quarto.yml b/docs/_quarto.yml
index e0dda9d6..eb791b59 100644
--- a/docs/_quarto.yml
+++ b/docs/_quarto.yml
@@ -59,6 +59,7 @@ quartodoc:
       contents:
         - estimation.estimation.feols
         - estimation.estimation.fepois
+        - estimation.estimation.feglm
         - did.estimation.did2s
         - did.estimation.lpdid
         - did.estimation.event_study
@@ -71,6 +72,10 @@ quartodoc:
         - estimation.feols_.Feols
         - estimation.fepois_.Fepois
         - estimation.feiv_.Feiv
+        - estimation.feglm_.Feglm
+        - estimation.felogit_.Felogit
+        - estimation.feprobit_.Feprobit
+        - estimation.fegaussian_.Fegaussian
         - estimation.feols_compressed_.FeolsCompressed
         #- did.did.DID
         #- did.did2s.DID2s
diff --git a/docs/_sidebar.yml b/docs/_sidebar.yml
index 80653ddb..304bf5da 100644
--- a/docs/_sidebar.yml
+++ b/docs/_sidebar.yml
@@ -5,6 +5,7 @@ website:
     - contents:
       - reference/estimation.estimation.feols.qmd
       - reference/estimation.estimation.fepois.qmd
+      - reference/estimation.estimation.feglm.qmd
       - reference/did.estimation.did2s.qmd
       - reference/did.estimation.lpdid.qmd
       - reference/did.estimation.event_study.qmd
@@ -15,6 +16,10 @@ website:
       - reference/estimation.feols_.Feols.qmd
       - reference/estimation.fepois_.Fepois.qmd
       - reference/estimation.feiv_.Feiv.qmd
+      - reference/estimation.feglm_.Feglm.qmd
+      - reference/estimation.felogit_.Felogit.qmd
+      - reference/estimation.feprobit_.Feprobit.qmd
+      - reference/estimation.fegaussian_.Fegaussian.qmd
       - reference/estimation.feols_compressed_.FeolsCompressed.qmd
       section: Estimation Classes
     - contents:
diff --git a/pyfixest/estimation/estimation.py b/pyfixest/estimation/estimation.py
index 64485022..8e69eaa8 100644
--- a/pyfixest/estimation/estimation.py
+++ b/pyfixest/estimation/estimation.py
@@ -641,8 +641,8 @@ def fepois(
     fit = pf.fepois("Y ~ X1 + X2 | f1 + f2", data)
     fit.summary()
     ```
-    For more examples, please take a look at the documentation of the `feols()`
-    function.
+
+    For more examples on the use of other function arguments, please take a look at the documentation of the [feols()](https://py-econometrics.github.io/pyfixest/reference/estimation.estimation.feols.html#pyfixest.estimation.estimation.feols) function.
     """
     if separation_check is None:
         separation_check = ["fe"]
@@ -859,13 +859,40 @@ def feglm(
 
     ```{python}
     import pyfixest as pf
+    import numpy as np
 
-    data = pf.get_data(model = "Fepois")
-    fit = pf.fepois("Y ~ X1 + X2 | f1 + f2", data)
-    fit.summary()
+    data = pf.get_data()
+    data["Y"] = np.where(data["Y"] > 0, 1, 0)
+    data["f1"] = np.where(data["f1"] > data["f1"].median(), "group1", "group2")
+
+    fit_probit = pf.feglm("Y ~ X1*f1", data, family = "probit")
+    fit_logit = pf.feglm("Y ~ X1*f1", data, family = "logit")
+    fit_gaussian = pf.feglm("Y ~ X1*f1", data, family = "gaussian")
+
+    pf.etable([fit_probit, fit_logit, fit_gaussian])
+    ```
+
+    `PyFixest` integrates with the [marginaleffects](https://marginaleffects.com/bonus/python.html) package. For example, to compute average marginal effects
+    for the probit model above, you can use the following code:
+
+    ```{python}
+    # we load polars as marginaleffects outputs pl.DataFrame's
+    import polars as pl
+    from marginaleffects import avg_slopes
+    pl.concat([avg_slopes(model, variables  = "X1") for model in [fit_probit, fit_logit, fit_gaussian]])
+    ```
+
+    We can also compute marginal effects by group (group average marginal effects):
+
+    ```{python}
+    avg_slopes(fit_probit, variables  = "X1", by = "f1")
     ```
-    For more examples, please take a look at the documentation of the `feols()`
+
+    We find homogeneous effects by "f1" in the probit model.
+
+    For more examples of other function arguments, please take a look at the documentation of the [feols()](https://py-econometrics.github.io/pyfixest/reference/estimation.estimation.feols.html#pyfixest.estimation.estimation.feols)
     function.
+
     """
     if family not in ["logit", "probit", "gaussian"]:
         raise ValueError(