diff --git a/bambi/backend/links.py b/bambi/backend/links.py
index 9e5af25d0..f723a8aaa 100644
--- a/bambi/backend/links.py
+++ b/bambi/backend/links.py
@@ -3,7 +3,7 @@
 
 
 def probit(x):
-    """Probit function that ensures result is in (0, 1)"""
+    """Probit function that ensures result is in (0, 1)."""
     eps = np.finfo(float).eps
     result = 0.5 + 0.5 * pt.erf(x / pt.sqrt(2))
     result = pt.switch(pt.eq(result, 0), eps, result)
@@ -13,7 +13,7 @@ def probit(x):
 
 
 def cloglog(x):
-    """Cloglog function that ensures result is in (0, 1)"""
+    """Cloglog function that ensures result is in (0, 1)."""
     eps = np.finfo(float).eps
     result = 1 - pt.exp(-pt.exp(x))
     result = pt.switch(pt.eq(result, 0), eps, result)
@@ -23,7 +23,7 @@ def cloglog(x):
 
 
 def logit(x):
-    """Logit function that ensures result is in (0, 1)"""
+    """Logit function that ensures result is in (0, 1)."""
     eps = np.finfo(float).eps
     result = pt.sigmoid(x)
     result = pt.switch(pt.eq(result, 0), eps, result)
diff --git a/bambi/backend/model_components.py b/bambi/backend/model_components.py
index 686959fa5..09e237ccf 100644
--- a/bambi/backend/model_components.py
+++ b/bambi/backend/model_components.py
@@ -66,7 +66,7 @@ def build_intercept(self, bmb_model):
             self.output += InterceptTerm(self.component.intercept_term).build(bmb_model)
 
     def build_offsets(self):
-        """Add intercept term to the PyMC model.
+        """Add intercept term to the PyMC model
 
         We have linear predictors of the form 'X @ b + Z @ u'. This is technically part of
         'X @ b' but it is added separately for convenience reasons.
@@ -75,10 +75,10 @@ def build_offsets(self):
             self.output += offset.data.squeeze()
 
     def build_common_terms(self, pymc_backend, bmb_model):
-        """Add common (fixed) terms to the PyMC model.
+        """Add common (fixed) terms to the PyMC model
 
         We have linear predictors of the form 'X @ b + Z @ u'.
-        This creates the 'b' parameter vector in PyMC, computes `X @ b`, and adds it to ``self.mu``.
+        This creates the 'b' parameter vector in PyMC, computes 'X @ b', and adds it to ``self.mu``.
 
         Parameters
         ----------
@@ -115,7 +115,7 @@ def build_common_terms(self, pymc_backend, bmb_model):
             self.output += pt.dot(data, coefs)
 
     def build_hsgp_terms(self, pymc_backend, bmb_model):
-        """Add HSGP (Hilbert-Space Gaussian Process approximation) terms to the PyMC model.
+        """Add HSGP (Hilbert-Space Gaussian Process approximation) terms to the PyMC model
 
         The linear predictor 'X @ b + Z @ u' can be augmented with non-parametric HSGP terms
         'f(x)'. This creates the 'f(x)' and adds it ``self.output``.
@@ -128,10 +128,10 @@ def build_hsgp_terms(self, pymc_backend, bmb_model):
             self.output += hsgp_term.build(bmb_model)
 
     def build_group_specific_terms(self, pymc_backend, bmb_model):
-        """Add group-specific (random or varying) terms to the PyMC model.
+        """Add group-specific (random or varying) terms to the PyMC model
 
         We have linear predictors of the form 'X @ b + Z @ u'.
-        This creates the 'u' parameter vector in PyMC, computes `Z @ u`, and adds it to
+        This creates the 'u' parameter vector in PyMC, computes 'Z @ u', and adds it to
         ``self.output``.
         """
         for term in self.component.group_specific_terms.values():
diff --git a/bambi/backend/pymc.py b/bambi/backend/pymc.py
index 354f08ce3..9c46c045a 100644
--- a/bambi/backend/pymc.py
+++ b/bambi/backend/pymc.py
@@ -48,7 +48,7 @@ def __init__(self):
         self.components = {}
 
     def build(self, spec):
-        """Compile the PyMC model from an abstract model specification.
+        """Compile the PyMC model from an abstract model specification
 
         Parameters
         ----------
@@ -132,7 +132,7 @@ def build_response(self, spec):
         response_component.build_response(self, spec)
 
     def build_potentials(self, spec):
-        """Add potentials to the PyMC model.
+        """Add potentials to the PyMC model
 
         Potentials are arbitrary quantities that are added to the model log likelihood.
         See 'Factor Potentials' in
@@ -313,7 +313,7 @@ def _run_vi(self, **kwargs):
         return self.vi_approx
 
     def _run_laplace(self, draws, omit_offsets, include_mean):
-        """Fit a model using a Laplace approximation.
+        """Fit a model using a Laplace approximation
 
         Mainly for pedagogical use, provides reasonable results for approximately
         Gaussian posteriors. The approximation can be very poor for some models
@@ -369,12 +369,12 @@ def distributional_components(self):
 
 
 def _posterior_samples_to_idata(samples, model):
-    """Create InferenceData from samples.
+    """Create InferenceData from samples
 
     Parameters
     ----------
     samples: array
-        Posterior samples
+        Posterior samples.
     model: PyMC model
 
     Returns
diff --git a/bambi/backend/terms.py b/bambi/backend/terms.py
index 9f9dff549..031398b06 100644
--- a/bambi/backend/terms.py
+++ b/bambi/backend/terms.py
@@ -164,7 +164,7 @@ def name(self):
 
 
 class InterceptTerm:
-    """Representation of an intercept term in a PyMC model.
+    """Representation of an intercept term in a PyMC model
 
     Parameters
     ----------
@@ -194,7 +194,7 @@ def name(self):
 
 
 class ResponseTerm:
-    """Representation of a response term in a PyMC model.
+    """Representation of a response term in a PyMC model
 
     Parameters
     ----------
@@ -209,19 +209,19 @@ def __init__(self, term, family):
         self.family = family
 
     def build(self, pymc_backend, bmb_model):
-        """Create and return the response distribution for the PyMC model.
+        """Create and return the response distribution for the PyMC model
 
         Parameters
         ----------
         pymc_backend : bambi.backend.PyMCModel
             The object with all the backend information
         bmb_model : bambi.Model
-            The Bambi model instance
+            The Bambi model instance.
 
         Returns
         -------
         dist : pm.Distribution
-            The response distribution
+            The response distribution.
         """
         data = np.squeeze(self.term.data)
         parent = self.family.likelihood.parent
@@ -504,7 +504,7 @@ def get_covariance_functions(self):
         Returns
         -------
         Sequence[pm.gp.Covariance]
-            A covariance function that can be used with a GP in PyMC
+            A covariance function that can be used with a GP in PyMC.
         """
 
         # Get the callable that creates the function
diff --git a/bambi/backend/utils.py b/bambi/backend/utils.py
index 6827099d9..6792b1591 100644
--- a/bambi/backend/utils.py
+++ b/bambi/backend/utils.py
@@ -57,7 +57,7 @@ def get_linkinv(link, invlinks):
     Returns
     -------
         callable
-        The link function
+        The link function.
     """
     # If the name is in the backend, get it from there
     if link.name in invlinks:
@@ -97,7 +97,7 @@ def make_weighted_logp(dist: pm.Distribution):
 
     Returns
     -------
-    A function that computes the weighted logp
+    A function that computes the weighted logp.
     """
 
     def logp(value, *dist_params, weights):
@@ -108,7 +108,7 @@ def logp(value, *dist_params, weights):
 
 
 def get_dist_args(dist: pm.Distribution) -> list[str]:
-    """Get the argument names of a PyMC distribution.
+    """Get the argument names of a PyMC distribution
 
     The argument names are the names of the parameters of the distribution.
 
diff --git a/bambi/families/family.py b/bambi/families/family.py
index 9f1375ced..5466215a5 100644
--- a/bambi/families/family.py
+++ b/bambi/families/family.py
@@ -9,7 +9,7 @@
 
 
 class Family:
-    """A specification of model family.
+    """A specification of model family
 
     Parameters
     ----------
@@ -110,7 +110,7 @@ def posterior_predictive(self, model, posterior, **kwargs):
 
         This function works for almost all the families. It grabs the draws for the parameters
         needed in the response distribution, and then gets samples from the posterior predictive
-        distribution using `pm.draw()`. It won't work when the response distribution requires
+        distribution using ``pm.draw()``. It won't work when the response distribution requires
         parameters that are not available in `posterior`.
 
         Parameters
@@ -129,7 +129,7 @@ def posterior_predictive(self, model, posterior, **kwargs):
         Returns
         -------
         xr.DataArray
-            A data array with the draws from the posterior predictive distribution
+            A data array with the draws from the posterior predictive distribution.
         """
         response_dist = get_response_dist(model.family)
         params = model.family.likelihood.params
@@ -213,12 +213,12 @@ def get_response_dist(family):
     Parameters
     ----------
     family : bambi.Family
-        The family for which the response distribution is wanted
+        The family for which the response distribution is wanted.
 
     Returns
     -------
     pm.Distribution
-        The response distribution
+        The response distribution.
     """
     mapping = {"Cumulative": pm.Categorical, "StoppingRatio": pm.Categorical}
 
@@ -234,25 +234,25 @@ def get_response_dist(family):
 def expand_array(x, ndim):
     """Add dimensions to an array to match the number of desired dimensions
 
-    If x.ndim < ndim, it adds ndim - x.ndim dimensions after the last axis. If not, it is left
+    If ``x.ndim < ndim``, it adds ``ndim - x.ndim`` dimensions after the last axis. If not, it is left
     untouched.
 
-    For example, if we have a normal regression model with n = 1000, chains = 2, and draws = 500
-    the shape of the draws of mu will be (2, 500, 1000) but the shape of the draws of sigma will be
-    (2, 500). This function makes sure the shape of the draws of sigma is (2, 500, 1) which is
-    comaptible with (2, 500, 1000).
+    For example, if we have a normal regression model with ``n = 1000``, ``chains = 2``, and
+    ``draws = 500`` the shape of the draws of mu will be ``(2, 500, 1000)`` but the shape of the
+    draws of sigma will be ``(2, 500)``. This function makes sure the shape of the draws of
+    sigma is ``(2, 500, 1)`` which is comaptible with ``(2, 500, 1000)``.
 
     Parameters
     ----------
     x : np.ndarray
         The array
     ndim : int
-        The number of desired dimensions
+        The number of desired dimensions.
 
     Returns
     -------
     np.ndarray
-        The array with the expanded dimensions
+        The array with the expanded dimensions.
     """
     if x.ndim == ndim:
         return x
diff --git a/bambi/families/likelihood.py b/bambi/families/likelihood.py
index 99c8af311..92aebaaa5 100644
--- a/bambi/families/likelihood.py
+++ b/bambi/families/likelihood.py
@@ -29,7 +29,7 @@
 
 
 class Likelihood:
-    """Representation of a Likelihood function for a Bambi model.
+    """Representation of a Likelihood function for a Bambi model
 
     Notes:
     * ``parent`` must be in ``params``
diff --git a/bambi/families/link.py b/bambi/families/link.py
index a3abd8fee..0b194eab9 100644
--- a/bambi/families/link.py
+++ b/bambi/families/link.py
@@ -8,7 +8,7 @@
 
 
 def force_within_unit_interval(x):
-    """Make sure data in unit interval is in (0, 1)"""
+    """Make sure data in unit interval is in (0, 1)."""
     eps = np.finfo(float).eps
     x[x == 0] = eps
     x[x == 1] = 1 - eps
@@ -16,7 +16,7 @@ def force_within_unit_interval(x):
 
 
 def force_greater_than_zero(x):
-    """Make sure data in positive reals is in (0, infty)"""
+    """Make sure data in positive reals is in (0, infty)."""
     eps = np.finfo(float).eps
     x[x == 0] = eps
     return x
@@ -33,32 +33,32 @@ def cloglog(mu):
 
 
 def invcloglog(eta):
-    """Inverse of the cloglog function that ensures result is in (0, 1)"""
+    """Inverse of the cloglog function that ensures result is in (0, 1)."""
     result = 1 - np.exp(-np.exp(eta))
     return force_within_unit_interval(result)
 
 
 def probit(mu):
-    """Probit function that ensures the input is in (0, 1)"""
+    """Probit function that ensures the input is in (0, 1)."""
     mu = force_within_unit_interval(mu)
     return 2**0.5 * special.erfinv(2 * mu - 1)  # pylint: disable=no-member
 
 
 def invprobit(eta):
-    """Inverse of the probit function that ensures result is in (0, 1)"""
+    """Inverse of the probit function that ensures result is in (0, 1)."""
     result = 0.5 + 0.5 * special.erf(eta / 2**0.5)  # pylint: disable=no-member
     return force_within_unit_interval(result)
 
 
 def expit(eta):
-    """Expit function that ensures result is in (0, 1)"""
+    """Expit function that ensures result is in (0, 1)."""
     result = special.expit(eta)  # pylint: disable=no-member
     result = force_within_unit_interval(result)
     return result
 
 
 def logit(mu):
-    """Logit function that ensures the input is in (0, 1)"""
+    """Logit function that ensures the input is in (0, 1)."""
     mu = force_within_unit_interval(mu)
     return special.logit(mu)  # pylint: disable=no-member
 
@@ -115,7 +115,7 @@ def link_not_implemented(*args, **kwargs):
 
 
 class Link:
-    """Representation of a link function.
+    """Representation of a link function
 
     This object contains two main functions. One is the link function itself, the function
     that maps values in the response scale to the linear predictor, and the other is the inverse
diff --git a/bambi/families/univariate.py b/bambi/families/univariate.py
index 806e10d69..e311c91f9 100644
--- a/bambi/families/univariate.py
+++ b/bambi/families/univariate.py
@@ -55,7 +55,7 @@ def get_success_level(self, response):
 
 
 class Beta(UnivariateFamily):
-    """Beta Family
+    """Beta family
 
     It uses the mean (mu) and sample size (kappa) parametrization of the Beta distribution.
     """
@@ -465,9 +465,9 @@ class ZeroInflatedPoisson(UnivariateFamily):
 
 # pylint: disable = protected-access
 def get_success_level(term):
-    """Returns the success level of a categorical term.
+    """Returns the success level of a categorical term
 
-    Whenever the concept of "success level" does not apply, it returns None.
+    Whenever the concept of "success level" does not apply, it returns ``None``.
     """
     if term.kind != "categoric":
         return None
@@ -485,9 +485,9 @@ def get_success_level(term):
 
 # pylint: disable = protected-access
 def get_reference_level(term):
-    """Returns the reference level of a categorical term.
+    """Returns the reference level of a categorical term
 
-    Whenever the concept of "reference level" does not apply, it returns None.
+    Whenever the concept of "reference level" does not apply, it returns ``None``.
     """
     if term.kind != "categoric":
         return None
diff --git a/bambi/formula.py b/bambi/formula.py
index ebcd94cfe..aaeace93b 100644
--- a/bambi/formula.py
+++ b/bambi/formula.py
@@ -17,7 +17,7 @@ class Formula:
     formula : str
         A model description written using the formula syntax from the ``formulae`` library.
     *additionals : str
-        Additional formulas that describe the
+        Additional formulas that describe model parameters rather than a response variable.
     """
 
     def __init__(self, formula: str, *additionals: str):
@@ -31,7 +31,7 @@ def check_additionals(self, additionals: Sequence[str]):
         Parameters
         ----------
         additionals : Sequence[str]
-            Model formulas that describe model parameters rather than a response variable
+            Model formulas that describe model parameters rather than a response variable.
 
         Returns
         -------
@@ -53,9 +53,9 @@ def check_additional(self, additional: str):
         Raises
         ------
         ValueError
-            If the formula does not contain a response term
+            If the formula does not contain a response term.
         ValueError
-            If the response term is not a plain name
+            If the response term is not a plain name.
         """
         response = fm.model_description(additional).response
 
@@ -75,7 +75,7 @@ def get_all_formulas(self):
         Returns
         -------
         list
-            All the formulas in the instance
+            All the formulas in the instance.
         """
         return [self.main] + list(self.additionals)
 
@@ -97,13 +97,13 @@ def formula_has_intercept(formula: str) -> bool:
 
 
 def check_ordinal_formula(formula: Formula) -> Formula:
-    """Check if a supplied formula can be used with an ordinal model.
+    """Check if a supplied formula can be used with an ordinal model
 
     Ordinal models have the following constrains (for the moment):
-    * A single formula must be passed. This is because Bambi does not support modeling the
-    thresholds as a function of predictors.
-    * The intercept is omitted. This is to avoid non-identifiability issues between the intercept
-    and the thresholds.
+    * A single formula must be passed. This is because Bambi does not support
+    modeling the thresholds as a function of predictors.
+    * The intercept is omitted. This is to avoid non-identifiability issues
+    between the intercept and the thresholds.
     """
     if len(formula.additionals) > 0:
         raise ValueError("Ordinal families don't accept multiple formulas")
diff --git a/bambi/interpret/create_data.py b/bambi/interpret/create_data.py
index 9711b3200..26bf1c3bd 100644
--- a/bambi/interpret/create_data.py
+++ b/bambi/interpret/create_data.py
@@ -18,7 +18,7 @@
 
 def _pairwise_grid(data_dict: dict) -> pd.DataFrame:
     """Creates a pairwise grid (cartesian product) of data by using the
-    key-values of the dictionary.
+    key-values of the dictionary
 
     Parameters
     ----------
@@ -44,9 +44,11 @@ def _grid_level(
     kind: str,
 ) -> pd.DataFrame:
     """Creates a "grid" of data by using the covariates passed into the
-    `conditional` argument. Values for the grid are either: (1) computed
-    using a equally spaced grid, mean, and or mode (depending on the
-    covariate dtype), and (2) a user specified value or range of values.
+    `conditional` argument
+
+    Values for the grid are either: (1) computed using a equally spaced
+    grid, mean, and or mode (depending on the covariate dtype), and (2)
+    a user specified value or range of values.
 
     Parameters
     ----------
@@ -131,8 +133,9 @@ def _grid_level(
 
 
 def _differences_unit_level(variable_info: VariableInfo, kind: str) -> pd.DataFrame:
-    """Creates the data for unit-level contrasts by using the observed (empirical)
-    data. All covariates in the model are included in the data, except for the
+    """Creates the data for unit-level contrasts by using the observed (empirical) data
+
+    All covariates in the model are included in the data, except for the
     contrast predictor. The contrast predictor is replaced with either: (1) the
     default contrast value, or (2) the user specified contrast value.
 
@@ -204,7 +207,7 @@ def create_differences_data(
 
 def create_predictions_data(condition_info: ConditionalInfo, user_passed: bool) -> pd.DataFrame:
     """Creates either unit level or grid level data for 'predictions' depending
-    if the user passed covariates.
+    if the user passed covariates
 
     Parameters
     ----------
diff --git a/bambi/interpret/effects.py b/bambi/interpret/effects.py
index 7537b49af..f7a557941 100644
--- a/bambi/interpret/effects.py
+++ b/bambi/interpret/effects.py
@@ -33,7 +33,7 @@
 class ResponseInfo:
     """Stores metadata about the response variable for indexing data in az.InferenceData,
     computing uncertainty intervals, and creating the summary dataframe in
-    'PredictiveDifferences'.
+    'PredictiveDifferences'
 
     Parameters
     ----------
@@ -75,8 +75,9 @@ def __post_init__(self):
 @dataclass
 class Estimate:
     """Stores the mean and bounds (uncertainty interval) of 'comparisons' and
-    'slopes' estimates. Used in 'PredictiveDifferences' to store typed data
-    for the summary dataframe.
+    'slopes' estimates
+
+    Used in 'PredictiveDifferences' to store typed data for the summary dataframe.
 
     Parameters
     ----------
@@ -128,7 +129,7 @@ def __post_init__(self):
 class PredictiveDifferences:
     """Computes predictive differences and their uncertainty intervals for
     'comparisons' and 'slopes' effects and returns a summary dataframe of the
-    results.
+    results
 
     Parameters
     ----------
@@ -235,10 +236,11 @@ def get_estimate(
         prob: float = 0.94,
     ):
         """Obtain the effect ('comparisons' or 'slopes') estimate and uncertainty
-        interval using the posterior samples. First, the posterior samples are
-        subsetted by the variable of interest's values. Then, the effect is
-        computed for each pairwise combination of the variable of interest's
-        values.
+        interval using the posterior samples
+
+        First, the posterior samples are subsetted by the variable of interest's
+        values. Then, the effect is computed for each pairwise combination of the
+        variable of interest's values.
 
         Parameters
         ----------
@@ -323,7 +325,8 @@ def get_estimate(
 
     def get_summary_df(self, response_dim: np.ndarray) -> pd.DataFrame:
         """
-        Builds the summary dataframe for 'comparisons' and 'slopes' effects.
+        Builds the summary dataframe for 'comparisons' and 'slopes' effects
+
         There are four scenarios to consider:
 
             1.) If the effect kind is 'comparisons' and more than 2 values are being
@@ -399,7 +402,7 @@ def get_summary_df(self, response_dim: np.ndarray) -> pd.DataFrame:
 
     def average_by(self, variable: Union[bool, str]) -> pd.DataFrame:
         """Uses the original 'summary_df' to perform a marginal (if 'variable=True')
-        or group by average if covariate(s) are passed.
+        or group by average if covariate(s) are passed
 
         Parameters
         ----------
diff --git a/bambi/interpret/logs.py b/bambi/interpret/logs.py
index 13ddb4116..919a6d324 100644
--- a/bambi/interpret/logs.py
+++ b/bambi/interpret/logs.py
@@ -5,7 +5,7 @@
 
 def log_interpret_defaults(func):
     """
-    Decorator for functions that compute default values.
+    Decorator for functions that compute default values
 
     Logs output to console if 'bmb.config["INTERPRET_VERBOSE"] = True' and when
     default values are computed for the variable of interest, i.e., 'contrast'
diff --git a/bambi/interpret/plot_types.py b/bambi/interpret/plot_types.py
index 4d56e1f92..d509cc9c4 100644
--- a/bambi/interpret/plot_types.py
+++ b/bambi/interpret/plot_types.py
@@ -13,7 +13,7 @@ def plot_numeric(
     legend: bool = True,
     axes=None,
 ):
-    """Plotting of numeric data types.
+    """Plotting of numeric data types
 
     Parameters
     ----------
@@ -24,11 +24,11 @@ def plot_numeric(
         function.
     transforms : dict
         Transformations that are applied to each of the variables being plotted. The keys are the
-        name of the variables, and the values are functions to be applied. Defaults to `None`.
+        name of the variables, and the values are functions to be applied. Defaults to ``None``.
     legend : bool, optional
-        Whether to include a legend in the plot. Default to `True`.
+        Whether to include a legend in the plot. Default to ``True``.
     axes : np.ndarray, optional
-        Array of axes. Defaults to `None`.
+        Array of axes. Defaults to ``None``.
 
     Returns
     -------
@@ -158,7 +158,7 @@ def plot_numeric(
 
 
 def plot_categoric(covariates: Covariates, plot_data: pd.DataFrame, legend: bool = True, axes=None):
-    """Plotting of categorical data types.
+    """Plotting of categorical data types
 
     Parameters
     ----------
@@ -168,9 +168,9 @@ def plot_categoric(covariates: Covariates, plot_data: pd.DataFrame, legend: bool
         The data created by the `create_cap_data` or `create_comparisons_data`
         function.
     legend : bool, optional
-        Whether to include a legend in the plot. Default to `True`.
+        Whether to include a legend in the plot. Default to ``True``.
     axes : np.ndarray, optional
-        Array of axes. Defaults to `None`.
+        Array of axes. Defaults to ``None``.
 
     Returns
     -------
diff --git a/bambi/interpret/utils.py b/bambi/interpret/utils.py
index 1ecd057f6..c905c0f09 100644
--- a/bambi/interpret/utils.py
+++ b/bambi/interpret/utils.py
@@ -20,10 +20,11 @@
 class VariableInfo:
     """
     Stores information about the variable (covariate) passed into the 'contrast'
-    or 'wrt' argument for comparisons and slopes. Depending on the effect type
-    ('slopes' or 'comparisons'), the values attribute is either: (1) the values
-    passed with the 'contrast' or 'wrt' argument, or (2) default are values
-    computed by calling 'set_default_variable_values()'.
+    or 'wrt' argument for comparisons and slopes
+
+    Depending on the effect type ('slopes' or 'comparisons'), the values attribute
+    is either: (1) the values passed with the 'contrast' or 'wrt' argument, or (2)
+    default are values computed by calling 'set_default_variable_values()'.
 
     'VariableInfo' is used to create 'slopes' and 'comparisons' data as well as
     computing the estimates for the 'slopes' and 'comparisons' effects.
@@ -91,17 +92,18 @@ def epsilon_difference(self, x, eps) -> np.ndarray:
     def set_default_variable_values(self) -> np.ndarray:
         """
         Returns default values for the variable of interest ('contrast' and 'wrt')
-        for the 'slopes' and 'comparisons' effects depending on the dtype of the
-        variable of interest, effect type, and if self.grid is True. The scenarios
-        are described below:
+        for the 'slopes' and 'comparisons' effects.
+
+        Depends on the dtype of the variable of interest, effect type, and if
+        self.grid is ``True``. The scenarios are described below:
 
         If numeric dtype and kind is 'comparisons', the returned value is a
         centered difference based on the mean of `variable'.
 
         If numeric dtype and kind is 'slopes', the returned value is an epsilon
-        difference based on the mean of `variable'.
+        difference based on the mean of 'variable'.
 
-        If categoric dtype the returned value is the unique levels of `variable'.
+        If categoric dtype the returned value is the unique levels of 'variable'.
         """
         terms = get_model_terms(self.model)
         # get default values for each variable in the model
@@ -136,7 +138,7 @@ def set_default_variable_values(self) -> np.ndarray:
 class ConditionalInfo:
     """
     Stores information about the conditional (covariates) passed into the
-    'conditional' argument for 'comparisons' and 'slopes' effects.
+    'conditional' argument for 'comparisons' and 'slopes' effects
 
     'ConditionalInfo' is used to create 'slopes' and 'comparisons' data as well
     as computing the estimates for the 'slopes' and 'comparisons' effects.
diff --git a/bambi/model_components.py b/bambi/model_components.py
index f4691e5e2..b757ff3e5 100644
--- a/bambi/model_components.py
+++ b/bambi/model_components.py
@@ -13,7 +13,7 @@ class ConstantComponent:
     """Constant model components
 
     This is a component for a target parameter that has no predictors. This could be seen as
-    an intercept-only model for that parameter. For example, tis is the case for sigma when
+    an intercept-only model for that parameter. For example, this is the case for sigma when
     a non-distributional Normal linear regression model is used.
 
     Parameters
@@ -21,9 +21,9 @@ class ConstantComponent:
     name : str
         The name of the component. For example "sigma", "alpha", or "kappa".
     priors : bambi.priors.Prior
-        The prior distribution for the parameter
+        The prior distribution for the parameter.
     spec : bambi.Model
-        The Bambi model
+        The Bambi model.
     """
 
     def __init__(self, name, prior, spec):
@@ -55,7 +55,7 @@ class DistributionalComponent:
         auxiliary parameter (``"param"``). When ``"data"`` this is actually modeling the "parent"
         parameter of the family.
     spec : bambi.Model
-        The Bambi model
+        The Bambi model.
     """
 
     def __init__(self, design, priors, response_name, response_kind, spec):
@@ -555,5 +555,5 @@ def prepare_prior(prior, kind, auto_scale):
     elif isinstance(prior, Prior):
         prior.auto_scale = False
     else:
-        raise ValueError("'prior' must be instance of Prior or None.")
+        raise ValueError("'prior' must be instance of Prior or ``None``.")
     return prior
diff --git a/bambi/models.py b/bambi/models.py
index 5d5eafe09..8d8f4ba67 100644
--- a/bambi/models.py
+++ b/bambi/models.py
@@ -37,7 +37,7 @@
 
 
 class Model:
-    """Specification of model class.
+    """Specification of model class
 
     Parameters
     ----------
@@ -96,7 +96,7 @@ class Model:
         to the intercept of the centered data.
     extra_namespace : dict, optional
         Additional user supplied variables with transformations or data to include in the
-        environment where the formula is evaluated. Defaults to `None`.
+        environment where the formula is evaluated. Defaults to ``None``.
     """
 
     # pylint: disable=too-many-instance-attributes
@@ -240,7 +240,7 @@ def fit(
         random_seed=None,
         **kwargs,
     ):
-        """Fit the model using PyMC.
+        """Fit the model using PyMC
 
         Parameters
         ----------
@@ -348,20 +348,20 @@ def fit(
         )
 
     def build(self):
-        """Set up the model for sampling/fitting.
+        """Set up the model for sampling/fitting
 
         Creates an instance of the underlying PyMC model and adds all the necessary terms to it.
 
         Returns
         -------
-        None
+        ``None``.
         """
         self.backend = PyMCModel()
         self.backend.build(self)
         self.built = True
 
     def set_priors(self, priors=None, common=None, group_specific=None):
-        """Set priors for one or more existing terms.
+        """Set priors for one or more existing terms
 
         Parameters
         ----------
@@ -375,7 +375,7 @@ def set_priors(self, priors=None, common=None, group_specific=None):
 
         Returns
         -------
-        None
+        ``None``.
         """
         kwargs = dict(zip(["priors", "common", "group_specific"], [priors, common, group_specific]))
         self._added_priors.update(kwargs)
@@ -442,7 +442,7 @@ def _set_priors(self, priors=None, common=None, group_specific=None):
                         component.update_priors(prior)
 
     def _set_family(self, family, link):
-        """Set the Family of the model.
+        """Set the Family of the model
 
         Parameters
         ----------
@@ -460,7 +460,7 @@ def _set_family(self, family, link):
 
         Returns
         -------
-        None
+        ``None``.
         """
 
         # If string, get builtin family
@@ -495,7 +495,7 @@ def set_alias(self, aliases):
 
         Returns
         -------
-        None
+        ``None``.
         """
         if not isinstance(aliases, dict):
             raise ValueError(f"'aliases' must be a dictionary, not a {type(aliases)}.")
@@ -663,7 +663,7 @@ def plot_priors(
 
         Returns
         -------
-        axes: matplotlib axes
+        axes: matplotlib axes.
         """
         self._check_built()
 
diff --git a/bambi/priors/prior.py b/bambi/priors/prior.py
index d4c2da420..803e8e180 100644
--- a/bambi/priors/prior.py
+++ b/bambi/priors/prior.py
@@ -3,7 +3,7 @@
 
 
 class Prior:
-    """Abstract specification of a term prior.
+    """Abstract specification of a term prior
 
     Parameters
     ----------
@@ -27,7 +27,7 @@ def __init__(self, name, auto_scale=True, dist=None, **kwargs):
         self.dist = dist
 
     def update(self, **kwargs):
-        """Update the arguments of the prior with additional arguments.
+        """Update the arguments of the prior with additional arguments
 
         Parameters
         ----------
diff --git a/bambi/terms/offset.py b/bambi/terms/offset.py
index 05902fc7c..ea7cebc97 100644
--- a/bambi/terms/offset.py
+++ b/bambi/terms/offset.py
@@ -2,7 +2,7 @@
 
 
 class OffsetTerm(CommonTerm):
-    """Representation of a single offset term.
+    """Representation of a single offset term
 
     Parameters
     ----------
diff --git a/bambi/terms/utils.py b/bambi/terms/utils.py
index b0cc5bc36..60781dad1 100644
--- a/bambi/terms/utils.py
+++ b/bambi/terms/utils.py
@@ -3,12 +3,12 @@
 
 
 def is_single_component(term) -> bool:
-    """Determines if formulae term contains a single component"""
+    """Determines if formulae term contains a single component."""
     return hasattr(term, "components") and len(term.components) == 1
 
 
 def is_call_component(component) -> bool:
-    """Determines if formulae component is the result of a function call"""
+    """Determines if formulae component is the result of a function call."""
     return isinstance(component, fm.terms.call.Call)
 
 
@@ -23,7 +23,7 @@ def is_call_of_kind(call, kind):
 
 
 def is_censored_response(term):
-    """Determines if a formulae term represents a censored response"""
+    """Determines if a formulae term represents a censored response."""
     if not is_single_component(term):
         return False
     component = term.components[0]  # get the first (and single) component
@@ -33,7 +33,7 @@ def is_censored_response(term):
 
 
 def is_truncated_response(term):
-    """Determines if a formulae term represents a truncated response"""
+    """Determines if a formulae term represents a truncated response."""
     if not is_single_component(term):
         return False
     component = term.components[0]  # get the first (and single) component
@@ -43,7 +43,7 @@ def is_truncated_response(term):
 
 
 def is_weighted_response(term):
-    """Determines if a formulae term represents a weighted response"""
+    """Determines if a formulae term represents a weighted response."""
     if not is_single_component(term):
         return False
     component = term.components[0]  # get the first (and single) component
diff --git a/bambi/transformations.py b/bambi/transformations.py
index 03a489a35..e7629eb97 100644
--- a/bambi/transformations.py
+++ b/bambi/transformations.py
@@ -5,7 +5,7 @@
 
 
 def c(*args):  # pylint: disable=invalid-name
-    """Concatenate columns into a 2D NumPy Array"""
+    """Concatenate columns into a 2D NumPy Array."""
     return np.column_stack(args)
 
 
@@ -21,7 +21,7 @@ def censored(*args):
     upper bound, only if it's interval censoring, and the third argument contains the censoring
     statuses.
 
-    Valid censoring statuses are
+    Valid censoring statuses are:
 
     * "left": left censoring
     * "none": no censoring
@@ -199,49 +199,49 @@ def __call__(
         Parameters
         ----------
         m : int, Sequence[int], ndarray
-            The number of basis vectors. See `HSGP.reconciliate_shape` to see how it is
+            The number of basis vectors. See ``HSGP.reconciliate_shape`? to see how it is
             broadcasted/recycled.
         L : float, Sequence[float], Sequence[Sequence[float]], ndarray, optional
-            The boundary of the variable space. See `HSGP.reconciliate_shape` to see how it is
-            broadcasted/recycled. Defaults to `None`.
+            The boundary of the variable space. See ``HSGP.reconciliate_shape`` to see how it is
+            broadcasted/recycled. Defaults to ``None``.
         c : float, Sequence[float], Sequence[Sequence[float]], ndarray, optional
-            The proportion extension factor. Se `HSGP.reconciliate_shape` to see how it is
-            broadcasted/recycled. Defaults to `None`.
+            The proportion extension factor. Se ``HSGP.reconciliate_shape`` to see how it is
+            broadcasted/recycled. Defaults to ``None``.
         by : array-like, optional
             The values of a variable to group by. It is used to create a HSGP term by group.
-            Defaults to `None`.
+            Defaults to ``None``.
         cov : str, optional
             The name of the covariance function to use. Defaults to "ExpQuad".
         share_cov : bool, optional
-            Whether to share the same covariance function for every group. Defaults to `True`.
+            Whether to share the same covariance function for every group. Defaults to ``True``.
         scale : bool, optional
-            When `True`, the predictors are be rescaled such that the largest Euclidean
+            When ``True``, the predictors are be rescaled such that the largest Euclidean
             distance between two points is 1. This adjustment often improves the sampling speed and
             convergence. The rescaling also impacts the estimated length-scale parameters,
             which will resemble those of the scaled predictors rather than the original predictors
-            when `scale` is `True`. Defaults to `None`, which means the behavior depends on whether
-            custom priors are passed or not. If custom priors are used, `None` is translated to
-            `False`. If automatic priors are used, `None` is translated to `True`.
+            when `scale` is ``True``. Defaults to ``None``, which means the behavior depends on whether
+            custom priors are passed or not. If custom priors are used, ``None`` is translated to
+            `False`. If automatic priors are used, ``None`` is translated to ``True``.
         iso : bool, optional
             Determines whether to use an isotropic or non-isotropic Gaussian Process.
             If isotropic, the same level of smoothing is applied to all predictors,
             while non-isotropic GPs allow different levels of smoothing for individual predictors.
-            This parameter is ignored if only one predictor is supplied. Defaults to `True`.
+            This parameter is ignored if only one predictor is supplied. Defaults to ``True``.
         drop_first : bool, optional
-            Whether to ignore the first basis vector or not. Defaults to `False`.
+            Whether to ignore the first basis vector or not. Defaults to ``False``.
         centered : bool, optional
-            Whether to use the centered or the non-centered parametrization. Defaults to `False`.
+            Whether to use the centered or the non-centered parametrization. Defaults to ``False``.
 
         Returns
         -------
         values
             A NumPy array of shape (observations_n, variables_n) or
-            (observations_n, variables_n + 1) if `by` is not None.
+            (observations_n, variables_n + 1) if ``by`` is not ``None``.
 
         Raises
         ------
         ValueError
-            When both `L` and `c` are `None` or when both of them are not `None` at the same time.
+            When both ``L`` and ``c`` are ``None`` or when both of them are not ``None`` at the same time.
         """
         values = np.column_stack(x)
 
@@ -304,19 +304,19 @@ def recycle_parameter(value, variables_n: int, groups_n: int):
         """Reshapes a value considering the number of variables and groups
 
         Parameter values such as `m`, `L`, and `c` may be different for the different variables and
-        groups. Internally, the shape of these objects is always `(groups_n, variables_n)`.
+        groups. Internally, the shape of these objects is always ``(groups_n, variables_n)``.
         This method contains the logic used to map user supplied values, which may be of different
-        shape and nature, into an object of shape `(groups_n, variables_n)`.
+        shape and nature, into an object of shape ``(groups_n, variables_n)``.
 
         The behavior of the method depends on the type of `value` in the following way.
         If value is of type...
         * `int`: the same value is recycled for all variables and groups.
         * `Sequence[int]`: it represents the values by variable and it is recycled for all groups.
         * `Sequence[Sequence[int]]`: it represents the values by variable and by group and thus
-        no recycling applies. Must be of shape `(groups_n, variables_n)`.
-        * `ndarray`:
-            * If one dimensional, it behaves as `Sequence[int]`
-            * If two dimensional, it behaves as `Sequence[Sequence[int]]`
+        no recycling applies. Must be of shape ``(groups_n, variables_n)``.
+        * ``ndarray``:
+            * If one dimensional, it behaves as ``Sequence[int]``
+            * If two dimensional, it behaves as ``Sequence[Sequence[int]]``
         """
         value = np.asarray(value)
         shape = value.shape
@@ -339,17 +339,17 @@ def as_matrix(x):
     Parameters
     ----------
     x : np.ndarray
-        Array
+        Array.
 
     Returns
     -------
     np.ndarray
-        A two dimensional array
+        A two dimensional array.
 
     Raises
     ------
     ValueError
-        If the input has more than two dimensions
+        If the input has more than two dimensions.
     """
     x = np.atleast_1d(x)
     if x.ndim == 1:
@@ -367,7 +367,7 @@ def mean_by_group(values, group):
     values : np.ndarray
         A 2 dimensional array. Rows indicate observations and columns indicate different variables.
     group : sequence
-        A sequence that indicates to which group each observation belongs to. If `None`, then
+        A sequence that indicates to which group each observation belongs to. If ``None``, then
         no group exists.
 
     Returns
diff --git a/bambi/utils.py b/bambi/utils.py
index 4de92f208..a01c65f08 100644
--- a/bambi/utils.py
+++ b/bambi/utils.py
@@ -10,7 +10,7 @@
 
 
 def listify(obj):
-    """Wrap all non-list or tuple objects in a list.
+    """Wrap all non-list or tuple objects in a list
 
     Provides a simple way to accept flexible arguments.
     """
@@ -33,7 +33,7 @@ def multilinify(sequence: Sequence[str], sep: str = ",") -> str:
 
 
 def wrapify(string: str, width: int = 100, indentation: int = 2) -> str:
-    """Wraps long strings into multiple lines.
+    """Wraps long strings into multiple lines
 
     This function is used to print the model summary.
     """
@@ -49,9 +49,9 @@ def wrapify(string: str, width: int = 100, indentation: int = 2) -> str:
 
 
 def extract_argument_names(expr, accepted_funcs):
-    """Extract the names of the arguments passed to a function.
+    """Extract the names of the arguments passed to a function
 
-    This is used to extract the labels from function calls such as `c(y1, y2, y3, y3)`.
+    This is used to extract the labels from function calls such as ``c(y1, y2, y3, y3)``.
 
     Parameters
     ----------
@@ -63,7 +63,7 @@ def extract_argument_names(expr, accepted_funcs):
     Returns
     -------
     list
-        If all criteria are met, the names of the arguments. Otherwise it returns None.
+        If all criteria are met, the names of the arguments. Otherwise it returns ``None``.
     """
     # Extract the first thing in the body
     parsed_expr = ast.parse(expr).body[0]
@@ -93,7 +93,7 @@ def extract_argument_names(expr, accepted_funcs):
 
 
 def clean_formula_lhs(x):
-    """Remove the left hand side of a model formula and the tilde.
+    """Remove the left hand side of a model formula and the tilde
 
     Parameters
     ----------
@@ -103,7 +103,7 @@ def clean_formula_lhs(x):
     Returns
     -------
     str
-        The right hand side of the model formula
+        The right hand side of the model formula.
     """
     assert "~" in x
     position = x.find("~")
@@ -119,12 +119,12 @@ def get_auxiliary_parameters(family):
     Parameters
     ----------
     family : bambi.families.Family
-        The family
+        The family.
 
     Returns
     -------
     set
-        Names of auxiliary parameters in the family
+        Names of auxiliary parameters in the family.
     """
     return set(family.likelihood.params) - {family.likelihood.parent}
 
@@ -138,12 +138,12 @@ def get_aliased_name(term):
     Parameters
     ----------
     term : BaseTerm
-        The term
+        The term.
 
     Returns
     -------
     str
-        The aliased name
+        The aliased name.
     """
     if term.alias:
         return term.alias
@@ -151,17 +151,17 @@ def get_aliased_name(term):
 
 
 def is_single_component(term) -> bool:
-    """Determines if formulae term contains a single component"""
+    """Determines if formulae term contains a single component."""
     return hasattr(term, "components") and len(term.components) == 1
 
 
 def is_call_component(component) -> bool:
-    """Determines if formulae component is the result of a function call"""
+    """Determines if formulae component is the result of a function call."""
     return isinstance(component, fm.terms.call.Call)
 
 
 def is_stateful_transform(component):
-    """Determines if formulae call component is a stateful transformation"""
+    """Determines if formulae call component is a stateful transformation."""
     return component.call.stateful_transform is not None