From c3aebfecf975b415ce8259d26be0d2742aa42665 Mon Sep 17 00:00:00 2001 From: Quarto GHA Workflow Runner Date: Thu, 16 Nov 2023 22:13:09 +0000 Subject: [PATCH] Built site for gh-pages --- .nojekyll | 2 +- notebooks/ESCS_multiple_regression.html | 3 --- notebooks/Strack_RRR_re_analysis.html | 3 --- notebooks/alternative_links_binary.html | 3 --- notebooks/beta_regression.html | 3 --- notebooks/categorical_regression.html | 3 --- notebooks/circular_regression.html | 3 --- notebooks/distributional_models.html | 3 --- notebooks/hierarchical_binomial_bambi.html | 3 --- notebooks/hsgp_1d.html | 4 +--- notebooks/hsgp_2d.html | 5 +---- notebooks/logistic_regression.html | 3 --- notebooks/mister_p.html | 3 --- notebooks/model_comparison.html | 3 --- notebooks/multi-level_regression.html | 3 --- notebooks/negative_binomial.html | 3 --- notebooks/ordinal_regression.html | 3 --- notebooks/plot_comparisons.html | 3 --- notebooks/plot_predictions.html | 3 --- notebooks/plot_slopes.html | 3 --- notebooks/predict_new_groups.html | 3 --- notebooks/quantile_regression.html | 3 --- notebooks/radon_example.html | 3 --- notebooks/shooter_crossed_random_ANOVA.html | 3 --- notebooks/sleepstudy.html | 3 --- notebooks/splines_cherry_blossoms.html | 3 --- notebooks/survival_model.html | 3 --- notebooks/t-test.html | 3 --- notebooks/t_regression.html | 3 --- notebooks/wald_gamma_glm.html | 3 --- notebooks/zero_inflated_regression.html | 3 --- search.json | 4 ++-- 32 files changed, 5 insertions(+), 94 deletions(-) diff --git a/.nojekyll b/.nojekyll index 722c3413c..299fcc02a 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -1bd86105 \ No newline at end of file +1fe8792d \ No newline at end of file diff --git a/notebooks/ESCS_multiple_regression.html b/notebooks/ESCS_multiple_regression.html index c9cf2cd73..279e19d30 100644 --- a/notebooks/ESCS_multiple_regression.html +++ b/notebooks/ESCS_multiple_regression.html @@ -385,9 +385,6 @@

Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/Strack_RRR_re_analysis.html b/notebooks/Strack_RRR_re_analysis.html index 9816fe9ac..fb173365e 100644 --- a/notebooks/Strack_RRR_re_analysis.html +++ b/notebooks/Strack_RRR_re_analysis.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/alternative_links_binary.html b/notebooks/alternative_links_binary.html index b315a94c9..92eb67f17 100644 --- a/notebooks/alternative_links_binary.html +++ b/notebooks/alternative_links_binary.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/beta_regression.html b/notebooks/beta_regression.html index 1719d8657..539535247 100644 --- a/notebooks/beta_regression.html +++ b/notebooks/beta_regression.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/categorical_regression.html b/notebooks/categorical_regression.html index 411b5346c..dc3ebd757 100644 --- a/notebooks/categorical_regression.html +++ b/notebooks/categorical_regression.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/circular_regression.html b/notebooks/circular_regression.html index 198c432a9..55978428c 100644 --- a/notebooks/circular_regression.html +++ b/notebooks/circular_regression.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/distributional_models.html b/notebooks/distributional_models.html index 4cfd74a2a..a9afa56e0 100644 --- a/notebooks/distributional_models.html +++ b/notebooks/distributional_models.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/hierarchical_binomial_bambi.html b/notebooks/hierarchical_binomial_bambi.html index 02bec7fa8..602bf290f 100644 --- a/notebooks/hierarchical_binomial_bambi.html +++ b/notebooks/hierarchical_binomial_bambi.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/hsgp_1d.html b/notebooks/hsgp_1d.html index 2196804eb..428f46d30 100644 --- a/notebooks/hsgp_1d.html +++ b/notebooks/hsgp_1d.html @@ -387,9 +387,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • @@ -420,6 +417,7 @@

    Gaussian Processes

    This article demonstrates the how to use Bambi with Gaussian Processes with 1 dimensional predictors. Bambi supports Gaussian Processes through the approximation known as Hilbert Space Gaussian Processes (HSGP).

    HSGP is a framework that falls under the class of low-rank approximations that are based on forming a basis function approximation with \(m\) basis functions, where \(m\) is usually much less smaller than \(n\), the number of observations.

    For references see Hilbert Space Methods for Reduced-Rank Gaussian Process Regression and Practical Hilbert Space Approximate Bayesian Gaussian Processes for Probabilistic Programming.

    +

    If you prefer a video format, have a look at Introduction to Hilbert Space GPs in PyMC given by Bill Engels.

    from formulae import design_matrices
 
diff --git a/notebooks/hsgp_2d.html b/notebooks/hsgp_2d.html
index 892be74c8..6a96d08eb 100644
--- a/notebooks/hsgp_2d.html
+++ b/notebooks/hsgp_2d.html
@@ -386,9 +386,6 @@ 
    
   Plot Slopes
    -
  • - notebooks/survival_models.ipynb -
    • @@ -427,7 +424,7 @@

    On this page

    Gaussian Processes in 2D

    This article demonstrates how to use Bambi with Gaussian Processes with 2 dimensional predictors. Bambi supports Gaussian Processes through the low-rank approximation known as Hilbert Space Gaussian Processes. For references see Hilbert Space Methods for Reduced-Rank Gaussian Process Regression and Practical Hilbert Space Approximate Bayesian Gaussian Processes for Probabilistic Programming.

    -

    For a demonstration of Gaussian Processes in 1D together with a more in depth explanation see To Do.

    +

    If you prefer a video format, have a look at Introduction to Hilbert Space GPs in PyMC given by Bill Engels.

    import arviz as az
 import bambi as bmb
diff --git a/notebooks/logistic_regression.html b/notebooks/logistic_regression.html
index 9492f8e5e..c3f061c4e 100644
--- a/notebooks/logistic_regression.html
+++ b/notebooks/logistic_regression.html
@@ -386,9 +386,6 @@ 
    
   Plot Slopes
    -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/mister_p.html b/notebooks/mister_p.html index 2ad70c0c1..8612c272a 100644 --- a/notebooks/mister_p.html +++ b/notebooks/mister_p.html @@ -386,9 +386,6 @@

    Plot Slopes
    -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/model_comparison.html b/notebooks/model_comparison.html index 7f5abc80f..1a5f96c91 100644 --- a/notebooks/model_comparison.html +++ b/notebooks/model_comparison.html @@ -386,9 +386,6 @@

    Plot Slopes
    -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/multi-level_regression.html b/notebooks/multi-level_regression.html index 436ef5fc1..d337771d2 100644 --- a/notebooks/multi-level_regression.html +++ b/notebooks/multi-level_regression.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/negative_binomial.html b/notebooks/negative_binomial.html index b831a9b6b..495b95415 100644 --- a/notebooks/negative_binomial.html +++ b/notebooks/negative_binomial.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/ordinal_regression.html b/notebooks/ordinal_regression.html index c322cc78a..26144dd22 100644 --- a/notebooks/ordinal_regression.html +++ b/notebooks/ordinal_regression.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/plot_comparisons.html b/notebooks/plot_comparisons.html index cac439570..46bcf1a7e 100644 --- a/notebooks/plot_comparisons.html +++ b/notebooks/plot_comparisons.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/plot_predictions.html b/notebooks/plot_predictions.html index 7adeca8bf..e9a1d4c4d 100644 --- a/notebooks/plot_predictions.html +++ b/notebooks/plot_predictions.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/plot_slopes.html b/notebooks/plot_slopes.html index e08cd3f0d..1e203d100 100644 --- a/notebooks/plot_slopes.html +++ b/notebooks/plot_slopes.html @@ -385,9 +385,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/predict_new_groups.html b/notebooks/predict_new_groups.html index 87d5c1568..f80716667 100644 --- a/notebooks/predict_new_groups.html +++ b/notebooks/predict_new_groups.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/quantile_regression.html b/notebooks/quantile_regression.html index b3b2aa5c4..d9154ec64 100644 --- a/notebooks/quantile_regression.html +++ b/notebooks/quantile_regression.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/radon_example.html b/notebooks/radon_example.html index e3c631272..732d7e10c 100644 --- a/notebooks/radon_example.html +++ b/notebooks/radon_example.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/shooter_crossed_random_ANOVA.html b/notebooks/shooter_crossed_random_ANOVA.html index cc397e0a9..dfcc79b7c 100644 --- a/notebooks/shooter_crossed_random_ANOVA.html +++ b/notebooks/shooter_crossed_random_ANOVA.html @@ -385,9 +385,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/sleepstudy.html b/notebooks/sleepstudy.html index 2be9b653f..2b8f47eba 100644 --- a/notebooks/sleepstudy.html +++ b/notebooks/sleepstudy.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/splines_cherry_blossoms.html b/notebooks/splines_cherry_blossoms.html index e534b3ffd..0cd0a9f98 100644 --- a/notebooks/splines_cherry_blossoms.html +++ b/notebooks/splines_cherry_blossoms.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/survival_model.html b/notebooks/survival_model.html index 367d7b0df..dbe3e9d90 100644 --- a/notebooks/survival_model.html +++ b/notebooks/survival_model.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/t-test.html b/notebooks/t-test.html index 80f9bd6bf..7bc026a8e 100644 --- a/notebooks/t-test.html +++ b/notebooks/t-test.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/t_regression.html b/notebooks/t_regression.html index fccd6be37..cf39f72e3 100644 --- a/notebooks/t_regression.html +++ b/notebooks/t_regression.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/wald_gamma_glm.html b/notebooks/wald_gamma_glm.html index 7cb30f16a..17e046a42 100644 --- a/notebooks/wald_gamma_glm.html +++ b/notebooks/wald_gamma_glm.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/notebooks/zero_inflated_regression.html b/notebooks/zero_inflated_regression.html index 3625a03ce..ae558f570 100644 --- a/notebooks/zero_inflated_regression.html +++ b/notebooks/zero_inflated_regression.html @@ -386,9 +386,6 @@

    Plot Slopes -
  • - notebooks/survival_models.ipynb -
    • diff --git a/search.json b/search.json index 7b42a4e84..e5020f811 100644 --- a/search.json +++ b/search.json @@ -263,14 +263,14 @@ "href": "notebooks/hsgp_1d.html", "title": "Bambi", "section": "", - "text": "This article demonstrates the how to use Bambi with Gaussian Processes with 1 dimensional predictors. Bambi supports Gaussian Processes through the approximation known as Hilbert Space Gaussian Processes (HSGP).\nHSGP is a framework that falls under the class of low-rank approximations that are based on forming a basis function approximation with \\(m\\) basis functions, where \\(m\\) is usually much less smaller than \\(n\\), the number of observations.\nFor references see Hilbert Space Methods for Reduced-Rank Gaussian Process Regression and Practical Hilbert Space Approximate Bayesian Gaussian Processes for Probabilistic Programming.\n\nfrom formulae import design_matrices\n\nimport arviz as az\nimport bambi as bmb\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\n\nfrom bambi.plots import plot_cap\nfrom matplotlib.lines import Line2D\n\n\n\nLet’s get started simulating some data from a smooth function. The goal is to fit a normal likelihood model where a Gaussian process term contributes to the mean.\n\nrng = np.random.default_rng(seed=121195)\n\nsize = 100\nx = np.linspace(0, 50, size)\nb = 0.1 * rng.normal(size=6)\nsigma = 0.15\n\ndm = design_matrices(\"0 + bs(x, df=6, intercept=True)\", pd.DataFrame({\"x\": x}))\nX = np.array(dm.common)\nf = 10 * X @ b\ny = f + rng.normal(size=size) * sigma\ndf = pd.DataFrame({\"x\": x, \"y\": y})\n\nfig, ax = plt.subplots(figsize=(9, 6))\nax.scatter(x, y, s=30, alpha=0.8);\nax.plot(x, f, color=\"black\");\n\n\n\n\nNow let’s simply create and fit the model. We use the hsgp to initialize a HSGP term in the model formula. Notice we pass the variable x and values for two other arguments m and c that we’ll cover later.\n\nmodel = bmb.Model(\"y ~ 0 + hsgp(x, m=10, c=2)\", df)\nmodel\n\n Formula: y ~ 0 + hsgp(x, m=10, c=2)\n Family: gaussian\n Link: mu = identity\n Observations: 100\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, m=10, c=2)\n cov: ExpQuad\n sigma ~ Exponential(lam: 1.0)\n ell ~ InverseGamma(alpha: 3.0, beta: 2.0)\n \n Auxiliary parameters\n y_sigma ~ HalfStudentT(nu: 4.0, sigma: 0.2745)\n\n\nIn the model description we can see the contribution of the HSGP term. It consists of two things: the name of the covariance kernel and the priors for its parameters. In this case, it’s an Exponentiated Quadratic covariance kernel with parameters sigma (amplitude) and ell (lengthscale). The prior for the amplitude is Exponential(1) and the prior for the lengthscale is InverseGamma(3, 2).\n\nidata = model.fit(inference_method=\"nuts_numpyro\", random_seed=121195)\nprint(idata.sample_stats[\"diverging\"].sum().to_numpy())\n\n/home/tomas/anaconda3/envs/bambi_hsgp/lib/python3.10/site-packages/pymc/sampling/jax.py:39: UserWarning: This module is experimental.\n warnings.warn(\"This module is experimental.\")\n\n\nCompiling...\n\n\nNo GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n\n\nCompilation time = 0:00:02.804363\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:04.776557\nTransforming variables...\nTransformation time = 0:00:00.521686\n527\n\n\n\naz.plot_trace(idata, backend_kwargs={\"layout\": \"constrained\"});\n\n\n\n\nThe fit is horrible. To fix that we can use better priors. But before doing that, it’s important to note that HSGP terms have a unique characteristic in that they do not receive priors themselves. Rather, the associated parameters of an HSGP term, such as sigma and ell, are the ones that are assigned priors. Therefore, we need to assign the HSGP term a dictionary of priors instead of a single prior.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=2), # amplitude\n \"ell\": bmb.Prior(\"InverseGamma\", mu=10, sigma=1) # lengthscale\n}\n\n# This is the dictionary we pass to Bambi\npriors = {\n \"hsgp(x, m=10, c=2)\": prior_hsgp,\n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=10)\n}\nmodel = bmb.Model(\"y ~ 0 + hsgp(x, m=10, c=2)\", df, priors=priors)\nmodel\n\n Formula: y ~ 0 + hsgp(x, m=10, c=2)\n Family: gaussian\n Link: mu = identity\n Observations: 100\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, m=10, c=2)\n cov: ExpQuad\n sigma ~ Exponential(lam: 2.0)\n ell ~ InverseGamma(mu: 10.0, sigma: 1.0)\n \n Auxiliary parameters\n y_sigma ~ HalfNormal(sigma: 10.0)\n\n\nNotice the priors were updated in the model summary. Now we’re ready to fit the model!\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9, random_seed=121195)\nprint(idata.sample_stats[\"diverging\"].sum().to_numpy())\n\nCompiling...\nCompilation time = 0:00:02.378503\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:05.336123\nTransforming variables...\nTransformation time = 0:00:00.174204\n7\n\n\n\naz.plot_trace(idata, backend_kwargs={\"layout\": \"constrained\"});\n\n\n\n\nThe marginal posteriors look somehow better, but we still have lots of divergences. What else can we do? Change the parametrization!\nThe hsgp() function has a centered argument which defaults to False and thus Bambi uses a non-centered parametrization by default. But we can change that actually. Let’s try it!\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=2), # amplitude\n \"ell\": bmb.Prior(\"InverseGamma\", mu=10, sigma=1) # lengthscale\n}\n\n# This is the dictionary we pass to Bambi\npriors = {\n \"hsgp(x, m=10, c=2, centered=True)\": prior_hsgp,\n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=10)\n}\nmodel = bmb.Model(\"y ~ 0 + hsgp(x, m=10, c=2, centered=True)\", df, priors=priors)\nmodel\n\n Formula: y ~ 0 + hsgp(x, m=10, c=2, centered=True)\n Family: gaussian\n Link: mu = identity\n Observations: 100\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, m=10, c=2, centered=True)\n cov: ExpQuad\n sigma ~ Exponential(lam: 2.0)\n ell ~ InverseGamma(mu: 10.0, sigma: 1.0)\n \n Auxiliary parameters\n y_sigma ~ HalfNormal(sigma: 10.0)\n\n\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9, random_seed=121195)\nprint(idata.sample_stats[\"diverging\"].sum().to_numpy())\n\nCompiling...\nCompilation time = 0:00:02.560797\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:04.839103\nTransforming variables...\nTransformation time = 0:00:00.028475\n0\n\n\n\naz.plot_trace(idata, backend_kwargs={\"layout\": \"constrained\"});\n\n\n\n\nAwesome! That looks much better now.\nWe still get all the nice things from Bambi when using GPs. An example of this is the plot_cap() function which enables us to generate a visualization of the adjusted mean with credible bands automatically.\n\nfig, ax = plt.subplots(figsize=(9, 6))\nax.scatter(df[\"x\"], df[\"y\"], s=30, color=\"0.5\", alpha=0.5)\nplot_cap(model, idata, \"x\", ax=ax);\nax.set(xlabel=\"Predictor\", ylabel=\"Observed\");\n\n\n\n\nAnd on top of that, it’s possible to get draws from the posterior predictive distribution and plot credible bands for it. All we need is the .predict() method from the model class.\n\nnew_data = pd.DataFrame({\"x\": np.linspace(0, 50, num=500)})\nmodel.predict(idata, kind=\"pps\", data=new_data)\npps = idata.posterior_predictive[\"y\"].to_numpy().reshape(4000, 500)\nqts = np.quantile(pps, q=(0.025, 0.975), axis=0)\n\nfig, ax = plt.subplots(figsize=(9, 6))\nax.fill_between(new_data[\"x\"], qts[0], qts[1], color=\"C0\", alpha=0.6)\nax.scatter(df[\"x\"], df[\"y\"], s=30, color=\"C1\", alpha=0.9)\nax.plot(x, f, color=\"black\", ls=\"--\")\nax.set(xlabel=\"Predictor\", ylabel=\"Observed\");\n\nhandles = [Line2D([], [], color=\"black\", ls=\"--\"), Line2D([], [], color=\"C0\")]\nlabels = [\"True curve\", \"Posterior predictive distribution\"]\nax.legend(handles, labels);\n\n\n\n\n\n\n\nhsgp() is a transformation that is available in the namespace where the model formula is evaluated. In plain english, hsgp() is like a function you can use in your model formulas. You don’t need to worry about the details, Bambi knows how to handle them.But if still you want to see the actual code, you can have a look at the implementation of the HSGP class in bambi/transformations.py.\nWhat users do need to care about is the arguments the hsgp() transformation support. There are a bunch of arguments that can be passed after the variable number of non-keyword arguments representing the variables of the HSGP contribution. Below is a brief overview of these arguments and their respective descriptions.\n\nm: The number of basis vectors\nL: The boundary of the variable space\nc: The proportion extension factor\nby: This argument specifies the values of a variable used for grouping. It is used to create a HSGP term by group. If left unspecified, the default value is None, which means that there is no group variable and all observations belong to the same group.\ncov: This argument specifies the name of the covariance function to be used. The default value is \"ExpQuad\".\nshare_cov: Determines whether the same covariance function is shared across all groups. This argument is relevant only when by is not None and the default value is True.\nscale: When set to 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.\niso: Determines whether to use an isotropic or non-isotropic Gaussian Process. With an isotropic GP, the same level of smoothing is applied to all predictors, while a anisotropic GP allows different levels of smoothing for individual predictors. Note that this argument is ignored if only one predictor is provided. The default value is True.\ndrop_first: Whether to exclude the first basis vector or not. The default value is False.\ncentered: Whether to use the centered or the non-centered parametrization. Defaults to False.\n\nThe parameters m, L and c are directly related to the basis vectors of the HSGP approximation. If you want to know more about m, L, and/or c, it’s recommended to have a look at the documentation of the HSGP class in PyMC.\n\nSo far, we showcased how to use m, c and centered. In the remainder of this article we’re going to see how by and share_cov are used when we add a GP contribution by groups.\n\n\n\nIn this section we fit a model with a HSGP contribution by levels of a categorical variable. The data was simulated with the gamSim() function from the R package {mgcv} by Simon Wood.\n\ndata = pd.read_csv(\"data/gam_data.csv\")\ndata[\"fac\"] = pd.Categorical(data[\"fac\"])\ndata.head()[[\"x2\", \"y\", \"fac\"]]\n\n\n\n\n\n \n \n \n x2\n y\n fac\n \n \n \n \n 0\n 0.497183\n 3.085274\n 3\n \n \n 1\n 0.196003\n -2.250410\n 2\n \n \n 2\n 0.958474\n 0.070548\n 3\n \n \n 3\n 0.972759\n -0.230454\n 1\n \n \n 4\n 0.755836\n 2.173497\n 2\n \n \n\n\n\n\nLet’s visualize x2 versus y for the different levels in fac.\n\nfig, ax = plt.subplots(figsize=(9, 5))\ncolors = [f\"C{i}\" for i in pd.Categorical(data[\"fac\"]).codes]\nax.scatter(data[\"x2\"], data[\"y\"], color=colors, alpha=0.6)\nax.set(xlabel=\"x2\", ylabel=\"y\");\n\n\n\n\nWe can observe the relation between x2 and y can be approximated by a smooth non-linear curve, for all groups.\nBelow, we create the model with Bambi. The biggest difference is that we’re passing by=fac in the hsgp() call. This is all we need to ask Bambi to create multiple GP contribution terms, one per group.\nAnother trick that was not shown previously is the usage of an alias. .set_alias() from the Model class allow us to have more readable and shorter names for the components of a model. As you’ll see below, it makes a huge difference when displaying summaries or visualizations for the parameters of the model.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"Exponential\", lam=3)\n}\npriors = {\n \"hsgp(x2, by=fac, m=12, c=1.5)\": prior_hsgp,\n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=1)\n}\nmodel = bmb.Model(\"y ~ 0 + hsgp(x2, by=fac, m=12, c=1.5)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x2, by=fac, m=12, c=1.5)\": \"hsgp\"})\nmodel\n\n Formula: y ~ 0 + hsgp(x2, by=fac, m=12, c=1.5)\n Family: gaussian\n Link: mu = identity\n Observations: 300\n Priors: \n target = mu\n HSGP contributions\n hsgp(x2, by=fac, m=12, c=1.5)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ Exponential(lam: 3.0)\n \n Auxiliary parameters\n y_sigma ~ HalfNormal(sigma: 1.0)\n\n\n\nmodel.build()\nmodel.graph()\n\n\n\n\nSee how nicer are the names for the HSGP contribution parameters with the alias!\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.95, random_seed=121195)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:03.565702\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:06.818602\nTransforming variables...\nTransformation time = 0:00:00.885410\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_weights\", \"hsgp_sigma\", \"hsgp_ell\", \"y_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nThis time we got no divergences and good mixing and nice convergence in our first try (or perhaps it wasn’t the first try!). One thing that stands out are the marginal posterior for some of the beta parameters (the weights of the basis). This may indicate our approximation is using more basis vectors than what’s really needed.\nNote: At this point we have used the term ‘basis vector’ several times. This concept is very close to the concept of ‘basis functions’. The difference is that the ‘basis vector’ is a ‘basis function’ already evaluated at a set of points. In this case, the set of points is made by the values of the numerical predictor x2.\nDo you remember how easy was it to use plot_cap() above? Should it be harder now? Well, the answer will surprise you: No!\nAll we need to do is passing a second variable name which is mapped to the color in the visualization. Voilà!\n\nfig, ax = plt.subplots(figsize = (9, 5))\ncolors = [f\"C{i}\" for i in pd.Categorical(data[\"fac\"]).codes]\nax.scatter(data[\"x2\"], data[\"y\"], color=colors, alpha=0.6)\nplot_cap(model, idata, [\"x2\", \"fac\"], ax=ax);\n\n\n\n\nWe can go one step further and modify the model to use different covariance functions for the different groups. For that purpose, we pass share_cov=False. As always, Bambi takes care of all the details.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=1),\n \"ell\": bmb.Prior(\"Exponential\", lam=3)\n}\npriors = {\n \"hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\": prior_hsgp,\n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=1)\n}\nmodel = bmb.Model(\"y ~ 0 + hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\": \"hsgp\"})\nmodel\n\n Formula: y ~ 0 + hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\n Family: gaussian\n Link: mu = identity\n Observations: 300\n Priors: \n target = mu\n HSGP contributions\n hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 1.0)\n ell ~ Exponential(lam: 3.0)\n \n Auxiliary parameters\n y_sigma ~ HalfNormal(sigma: 1.0)\n\n\n\nmodel.build()\nmodel.graph()\n\n\n\n\nHave a closer look at the model graph. See that the hsgp_sigma and hsgp_ell parameters are no longer scalar. There are three of each, one for each group.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.95, random_seed=121195)\n\nCompiling...\nCompilation time = 0:00:04.396845\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:07.743907\nTransforming variables...\nTransformation time = 0:00:00.519422\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"y_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nIn fact, we can see not all the groups have similar posteriors for the covariance function parameters when they are allowed to vary.\nBefore closing the article, it’s worth looking at a particular but not uncommon pattern when using the HSGP approximation. Let’s have a look at the posterior distributions for the weights of the basis.\n\naz.plot_trace(idata, var_names=[\"hsgp_weights\"], backend_kwargs={\"layout\": \"constrained\"});\n\n\n\n\nLooks like some distributions are extremely flat, and others are extremely tight around zero.\nTo investigate this further we can manually plot the posterior for the first J basis vectors and see what they look like.\n\nbasis_n = 6\nfig, axes = plt.subplots(3, 1, figsize = (7, 10))\nfor i in range(3):\n ax = axes[i]\n values = idata.posterior[\"hsgp_weights\"].sel({\"hsgp_by\": i + 1})\n for j in range(basis_n):\n az.plot_kde(\n values.sel({\"hsgp_weights_dim\": j}).to_numpy().flatten(), \n ax=ax, \n plot_kwargs={\"color\": f\"C{j}\"}\n );\n\n\n\n\nIndeed, we can see that, at least for the first group, the posterior of the weights start being too tight around zero when we consider the 6th basis vector. But the posteriors for the weights of the previous basis vectors look nice.\nTo confirm our thought, let’s increase the value of basis_n to 9 and see what happens.\n\nbasis_n = 9\nfig, axes = plt.subplots(3, 1, figsize = (7, 10))\nfor i in range(3):\n ax = axes[i]\n values = idata.posterior[\"hsgp_weights\"].sel({\"hsgp_by\": i + 1})\n for j in range(basis_n):\n az.plot_kde(\n values.sel({\"hsgp_weights_dim\": j}).to_numpy().flatten(), \n ax=ax, \n plot_kwargs={\"color\": f\"C{j}\"}\n );\n\n\n\n\nAlright, that’s too spiky! Nonetheless, we don’t see that happening for the third group yet, indicating the higher number of basis vectors is more appropriate for this group." + "text": "This article demonstrates the how to use Bambi with Gaussian Processes with 1 dimensional predictors. Bambi supports Gaussian Processes through the approximation known as Hilbert Space Gaussian Processes (HSGP).\nHSGP is a framework that falls under the class of low-rank approximations that are based on forming a basis function approximation with \\(m\\) basis functions, where \\(m\\) is usually much less smaller than \\(n\\), the number of observations.\nFor references see Hilbert Space Methods for Reduced-Rank Gaussian Process Regression and Practical Hilbert Space Approximate Bayesian Gaussian Processes for Probabilistic Programming.\nIf you prefer a video format, have a look at Introduction to Hilbert Space GPs in PyMC given by Bill Engels.\n\nfrom formulae import design_matrices\n\nimport arviz as az\nimport bambi as bmb\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\n\nfrom bambi.plots import plot_cap\nfrom matplotlib.lines import Line2D\n\n\n\nLet’s get started simulating some data from a smooth function. The goal is to fit a normal likelihood model where a Gaussian process term contributes to the mean.\n\nrng = np.random.default_rng(seed=121195)\n\nsize = 100\nx = np.linspace(0, 50, size)\nb = 0.1 * rng.normal(size=6)\nsigma = 0.15\n\ndm = design_matrices(\"0 + bs(x, df=6, intercept=True)\", pd.DataFrame({\"x\": x}))\nX = np.array(dm.common)\nf = 10 * X @ b\ny = f + rng.normal(size=size) * sigma\ndf = pd.DataFrame({\"x\": x, \"y\": y})\n\nfig, ax = plt.subplots(figsize=(9, 6))\nax.scatter(x, y, s=30, alpha=0.8);\nax.plot(x, f, color=\"black\");\n\n\n\n\nNow let’s simply create and fit the model. We use the hsgp to initialize a HSGP term in the model formula. Notice we pass the variable x and values for two other arguments m and c that we’ll cover later.\n\nmodel = bmb.Model(\"y ~ 0 + hsgp(x, m=10, c=2)\", df)\nmodel\n\n Formula: y ~ 0 + hsgp(x, m=10, c=2)\n Family: gaussian\n Link: mu = identity\n Observations: 100\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, m=10, c=2)\n cov: ExpQuad\n sigma ~ Exponential(lam: 1.0)\n ell ~ InverseGamma(alpha: 3.0, beta: 2.0)\n \n Auxiliary parameters\n y_sigma ~ HalfStudentT(nu: 4.0, sigma: 0.2745)\n\n\nIn the model description we can see the contribution of the HSGP term. It consists of two things: the name of the covariance kernel and the priors for its parameters. In this case, it’s an Exponentiated Quadratic covariance kernel with parameters sigma (amplitude) and ell (lengthscale). The prior for the amplitude is Exponential(1) and the prior for the lengthscale is InverseGamma(3, 2).\n\nidata = model.fit(inference_method=\"nuts_numpyro\", random_seed=121195)\nprint(idata.sample_stats[\"diverging\"].sum().to_numpy())\n\n/home/tomas/anaconda3/envs/bambi_hsgp/lib/python3.10/site-packages/pymc/sampling/jax.py:39: UserWarning: This module is experimental.\n warnings.warn(\"This module is experimental.\")\n\n\nCompiling...\n\n\nNo GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n\n\nCompilation time = 0:00:02.804363\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:04.776557\nTransforming variables...\nTransformation time = 0:00:00.521686\n527\n\n\n\naz.plot_trace(idata, backend_kwargs={\"layout\": \"constrained\"});\n\n\n\n\nThe fit is horrible. To fix that we can use better priors. But before doing that, it’s important to note that HSGP terms have a unique characteristic in that they do not receive priors themselves. Rather, the associated parameters of an HSGP term, such as sigma and ell, are the ones that are assigned priors. Therefore, we need to assign the HSGP term a dictionary of priors instead of a single prior.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=2), # amplitude\n \"ell\": bmb.Prior(\"InverseGamma\", mu=10, sigma=1) # lengthscale\n}\n\n# This is the dictionary we pass to Bambi\npriors = {\n \"hsgp(x, m=10, c=2)\": prior_hsgp,\n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=10)\n}\nmodel = bmb.Model(\"y ~ 0 + hsgp(x, m=10, c=2)\", df, priors=priors)\nmodel\n\n Formula: y ~ 0 + hsgp(x, m=10, c=2)\n Family: gaussian\n Link: mu = identity\n Observations: 100\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, m=10, c=2)\n cov: ExpQuad\n sigma ~ Exponential(lam: 2.0)\n ell ~ InverseGamma(mu: 10.0, sigma: 1.0)\n \n Auxiliary parameters\n y_sigma ~ HalfNormal(sigma: 10.0)\n\n\nNotice the priors were updated in the model summary. Now we’re ready to fit the model!\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9, random_seed=121195)\nprint(idata.sample_stats[\"diverging\"].sum().to_numpy())\n\nCompiling...\nCompilation time = 0:00:02.378503\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:05.336123\nTransforming variables...\nTransformation time = 0:00:00.174204\n7\n\n\n\naz.plot_trace(idata, backend_kwargs={\"layout\": \"constrained\"});\n\n\n\n\nThe marginal posteriors look somehow better, but we still have lots of divergences. What else can we do? Change the parametrization!\nThe hsgp() function has a centered argument which defaults to False and thus Bambi uses a non-centered parametrization by default. But we can change that actually. Let’s try it!\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=2), # amplitude\n \"ell\": bmb.Prior(\"InverseGamma\", mu=10, sigma=1) # lengthscale\n}\n\n# This is the dictionary we pass to Bambi\npriors = {\n \"hsgp(x, m=10, c=2, centered=True)\": prior_hsgp,\n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=10)\n}\nmodel = bmb.Model(\"y ~ 0 + hsgp(x, m=10, c=2, centered=True)\", df, priors=priors)\nmodel\n\n Formula: y ~ 0 + hsgp(x, m=10, c=2, centered=True)\n Family: gaussian\n Link: mu = identity\n Observations: 100\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, m=10, c=2, centered=True)\n cov: ExpQuad\n sigma ~ Exponential(lam: 2.0)\n ell ~ InverseGamma(mu: 10.0, sigma: 1.0)\n \n Auxiliary parameters\n y_sigma ~ HalfNormal(sigma: 10.0)\n\n\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9, random_seed=121195)\nprint(idata.sample_stats[\"diverging\"].sum().to_numpy())\n\nCompiling...\nCompilation time = 0:00:02.560797\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:04.839103\nTransforming variables...\nTransformation time = 0:00:00.028475\n0\n\n\n\naz.plot_trace(idata, backend_kwargs={\"layout\": \"constrained\"});\n\n\n\n\nAwesome! That looks much better now.\nWe still get all the nice things from Bambi when using GPs. An example of this is the plot_cap() function which enables us to generate a visualization of the adjusted mean with credible bands automatically.\n\nfig, ax = plt.subplots(figsize=(9, 6))\nax.scatter(df[\"x\"], df[\"y\"], s=30, color=\"0.5\", alpha=0.5)\nplot_cap(model, idata, \"x\", ax=ax);\nax.set(xlabel=\"Predictor\", ylabel=\"Observed\");\n\n\n\n\nAnd on top of that, it’s possible to get draws from the posterior predictive distribution and plot credible bands for it. All we need is the .predict() method from the model class.\n\nnew_data = pd.DataFrame({\"x\": np.linspace(0, 50, num=500)})\nmodel.predict(idata, kind=\"pps\", data=new_data)\npps = idata.posterior_predictive[\"y\"].to_numpy().reshape(4000, 500)\nqts = np.quantile(pps, q=(0.025, 0.975), axis=0)\n\nfig, ax = plt.subplots(figsize=(9, 6))\nax.fill_between(new_data[\"x\"], qts[0], qts[1], color=\"C0\", alpha=0.6)\nax.scatter(df[\"x\"], df[\"y\"], s=30, color=\"C1\", alpha=0.9)\nax.plot(x, f, color=\"black\", ls=\"--\")\nax.set(xlabel=\"Predictor\", ylabel=\"Observed\");\n\nhandles = [Line2D([], [], color=\"black\", ls=\"--\"), Line2D([], [], color=\"C0\")]\nlabels = [\"True curve\", \"Posterior predictive distribution\"]\nax.legend(handles, labels);\n\n\n\n\n\n\n\nhsgp() is a transformation that is available in the namespace where the model formula is evaluated. In plain english, hsgp() is like a function you can use in your model formulas. You don’t need to worry about the details, Bambi knows how to handle them.But if still you want to see the actual code, you can have a look at the implementation of the HSGP class in bambi/transformations.py.\nWhat users do need to care about is the arguments the hsgp() transformation support. There are a bunch of arguments that can be passed after the variable number of non-keyword arguments representing the variables of the HSGP contribution. Below is a brief overview of these arguments and their respective descriptions.\n\nm: The number of basis vectors\nL: The boundary of the variable space\nc: The proportion extension factor\nby: This argument specifies the values of a variable used for grouping. It is used to create a HSGP term by group. If left unspecified, the default value is None, which means that there is no group variable and all observations belong to the same group.\ncov: This argument specifies the name of the covariance function to be used. The default value is \"ExpQuad\".\nshare_cov: Determines whether the same covariance function is shared across all groups. This argument is relevant only when by is not None and the default value is True.\nscale: When set to 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.\niso: Determines whether to use an isotropic or non-isotropic Gaussian Process. With an isotropic GP, the same level of smoothing is applied to all predictors, while a anisotropic GP allows different levels of smoothing for individual predictors. Note that this argument is ignored if only one predictor is provided. The default value is True.\ndrop_first: Whether to exclude the first basis vector or not. The default value is False.\ncentered: Whether to use the centered or the non-centered parametrization. Defaults to False.\n\nThe parameters m, L and c are directly related to the basis vectors of the HSGP approximation. If you want to know more about m, L, and/or c, it’s recommended to have a look at the documentation of the HSGP class in PyMC.\n\nSo far, we showcased how to use m, c and centered. In the remainder of this article we’re going to see how by and share_cov are used when we add a GP contribution by groups.\n\n\n\nIn this section we fit a model with a HSGP contribution by levels of a categorical variable. The data was simulated with the gamSim() function from the R package {mgcv} by Simon Wood.\n\ndata = pd.read_csv(\"data/gam_data.csv\")\ndata[\"fac\"] = pd.Categorical(data[\"fac\"])\ndata.head()[[\"x2\", \"y\", \"fac\"]]\n\n\n\n\n\n \n \n \n x2\n y\n fac\n \n \n \n \n 0\n 0.497183\n 3.085274\n 3\n \n \n 1\n 0.196003\n -2.250410\n 2\n \n \n 2\n 0.958474\n 0.070548\n 3\n \n \n 3\n 0.972759\n -0.230454\n 1\n \n \n 4\n 0.755836\n 2.173497\n 2\n \n \n\n\n\n\nLet’s visualize x2 versus y for the different levels in fac.\n\nfig, ax = plt.subplots(figsize=(9, 5))\ncolors = [f\"C{i}\" for i in pd.Categorical(data[\"fac\"]).codes]\nax.scatter(data[\"x2\"], data[\"y\"], color=colors, alpha=0.6)\nax.set(xlabel=\"x2\", ylabel=\"y\");\n\n\n\n\nWe can observe the relation between x2 and y can be approximated by a smooth non-linear curve, for all groups.\nBelow, we create the model with Bambi. The biggest difference is that we’re passing by=fac in the hsgp() call. This is all we need to ask Bambi to create multiple GP contribution terms, one per group.\nAnother trick that was not shown previously is the usage of an alias. .set_alias() from the Model class allow us to have more readable and shorter names for the components of a model. As you’ll see below, it makes a huge difference when displaying summaries or visualizations for the parameters of the model.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"Exponential\", lam=3)\n}\npriors = {\n \"hsgp(x2, by=fac, m=12, c=1.5)\": prior_hsgp,\n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=1)\n}\nmodel = bmb.Model(\"y ~ 0 + hsgp(x2, by=fac, m=12, c=1.5)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x2, by=fac, m=12, c=1.5)\": \"hsgp\"})\nmodel\n\n Formula: y ~ 0 + hsgp(x2, by=fac, m=12, c=1.5)\n Family: gaussian\n Link: mu = identity\n Observations: 300\n Priors: \n target = mu\n HSGP contributions\n hsgp(x2, by=fac, m=12, c=1.5)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ Exponential(lam: 3.0)\n \n Auxiliary parameters\n y_sigma ~ HalfNormal(sigma: 1.0)\n\n\n\nmodel.build()\nmodel.graph()\n\n\n\n\nSee how nicer are the names for the HSGP contribution parameters with the alias!\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.95, random_seed=121195)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:03.565702\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:06.818602\nTransforming variables...\nTransformation time = 0:00:00.885410\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_weights\", \"hsgp_sigma\", \"hsgp_ell\", \"y_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nThis time we got no divergences and good mixing and nice convergence in our first try (or perhaps it wasn’t the first try!). One thing that stands out are the marginal posterior for some of the beta parameters (the weights of the basis). This may indicate our approximation is using more basis vectors than what’s really needed.\nNote: At this point we have used the term ‘basis vector’ several times. This concept is very close to the concept of ‘basis functions’. The difference is that the ‘basis vector’ is a ‘basis function’ already evaluated at a set of points. In this case, the set of points is made by the values of the numerical predictor x2.\nDo you remember how easy was it to use plot_cap() above? Should it be harder now? Well, the answer will surprise you: No!\nAll we need to do is passing a second variable name which is mapped to the color in the visualization. Voilà!\n\nfig, ax = plt.subplots(figsize = (9, 5))\ncolors = [f\"C{i}\" for i in pd.Categorical(data[\"fac\"]).codes]\nax.scatter(data[\"x2\"], data[\"y\"], color=colors, alpha=0.6)\nplot_cap(model, idata, [\"x2\", \"fac\"], ax=ax);\n\n\n\n\nWe can go one step further and modify the model to use different covariance functions for the different groups. For that purpose, we pass share_cov=False. As always, Bambi takes care of all the details.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=1),\n \"ell\": bmb.Prior(\"Exponential\", lam=3)\n}\npriors = {\n \"hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\": prior_hsgp,\n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=1)\n}\nmodel = bmb.Model(\"y ~ 0 + hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\": \"hsgp\"})\nmodel\n\n Formula: y ~ 0 + hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\n Family: gaussian\n Link: mu = identity\n Observations: 300\n Priors: \n target = mu\n HSGP contributions\n hsgp(x2, by=fac, m=12, c=1.5, share_cov=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 1.0)\n ell ~ Exponential(lam: 3.0)\n \n Auxiliary parameters\n y_sigma ~ HalfNormal(sigma: 1.0)\n\n\n\nmodel.build()\nmodel.graph()\n\n\n\n\nHave a closer look at the model graph. See that the hsgp_sigma and hsgp_ell parameters are no longer scalar. There are three of each, one for each group.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.95, random_seed=121195)\n\nCompiling...\nCompilation time = 0:00:04.396845\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:07.743907\nTransforming variables...\nTransformation time = 0:00:00.519422\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"y_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nIn fact, we can see not all the groups have similar posteriors for the covariance function parameters when they are allowed to vary.\nBefore closing the article, it’s worth looking at a particular but not uncommon pattern when using the HSGP approximation. Let’s have a look at the posterior distributions for the weights of the basis.\n\naz.plot_trace(idata, var_names=[\"hsgp_weights\"], backend_kwargs={\"layout\": \"constrained\"});\n\n\n\n\nLooks like some distributions are extremely flat, and others are extremely tight around zero.\nTo investigate this further we can manually plot the posterior for the first J basis vectors and see what they look like.\n\nbasis_n = 6\nfig, axes = plt.subplots(3, 1, figsize = (7, 10))\nfor i in range(3):\n ax = axes[i]\n values = idata.posterior[\"hsgp_weights\"].sel({\"hsgp_by\": i + 1})\n for j in range(basis_n):\n az.plot_kde(\n values.sel({\"hsgp_weights_dim\": j}).to_numpy().flatten(), \n ax=ax, \n plot_kwargs={\"color\": f\"C{j}\"}\n );\n\n\n\n\nIndeed, we can see that, at least for the first group, the posterior of the weights start being too tight around zero when we consider the 6th basis vector. But the posteriors for the weights of the previous basis vectors look nice.\nTo confirm our thought, let’s increase the value of basis_n to 9 and see what happens.\n\nbasis_n = 9\nfig, axes = plt.subplots(3, 1, figsize = (7, 10))\nfor i in range(3):\n ax = axes[i]\n values = idata.posterior[\"hsgp_weights\"].sel({\"hsgp_by\": i + 1})\n for j in range(basis_n):\n az.plot_kde(\n values.sel({\"hsgp_weights_dim\": j}).to_numpy().flatten(), \n ax=ax, \n plot_kwargs={\"color\": f\"C{j}\"}\n );\n\n\n\n\nAlright, that’s too spiky! Nonetheless, we don’t see that happening for the third group yet, indicating the higher number of basis vectors is more appropriate for this group." }, { "objectID": "notebooks/hsgp_2d.html", "href": "notebooks/hsgp_2d.html", "title": "Bambi", "section": "", - "text": "This article demonstrates how to use Bambi with Gaussian Processes with 2 dimensional predictors. Bambi supports Gaussian Processes through the low-rank approximation known as Hilbert Space Gaussian Processes. For references see Hilbert Space Methods for Reduced-Rank Gaussian Process Regression and Practical Hilbert Space Approximate Bayesian Gaussian Processes for Probabilistic Programming.\nFor a demonstration of Gaussian Processes in 1D together with a more in depth explanation see To Do.\n\nimport arviz as az\nimport bambi as bmb\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport pymc as pm\n\nThe goal of this notebook is to showcase Bambi’s support for Gaussian Processes on two-dimensional data using the HSGP approximation.\nTo achieve this, we begin by creating a matrix of coordinates that will serve as the locations where we measure the values of a continuous response variable.\n\nx1 = np.linspace(0, 10, 12)\nx2 = np.linspace(0, 10, 12)\nxx, yy = np.meshgrid(x1, x2)\nX = np.column_stack([xx.flatten(), yy.flatten()])\nX.shape\n\n(144, 2)\n\n\n\n\nIn modeling multi-dimensional data with a Gaussian Process, we must choose between using an isotropic or an anisotropic Gaussian Process. An isotropic GP applies the same degree of smoothing to all predictors and is rotationally invariant. On the other hand, an anisotropic GP assigns different degrees of smoothing to each predictor and is not rotationally invariant.\nFurthermore, as the hsgp() function allows for the creation of separate GP contribution terms for the levels of a categorical variable through its by argument, we also examine both single-group and multiple-group scenarios.\n\n\nWe create a covariance kernel using ExpQuad from the gp submodule in PyMC. Note that the lengthscale and amplitude for both dimensions are 2 and 1.2, respectively. Then, we simply use NumPy to get a random draw from the 144-dimensional multivariate normal distribution.\n\nrng = np.random.default_rng(1234)\n\nell = 2\ncov = 1.2 * pm.gp.cov.ExpQuad(2, ls=ell)\nK = cov(X).eval()\nmu = np.zeros(X.shape[0])\nprint(mu.shape, K.shape)\n\nf = rng.multivariate_normal(mu, K)\n\nfig, ax = plt.subplots()\nax.scatter(xx, yy, c=f, s=900, marker=\"s\");\n\n(144,) (144, 144)\n\n\n\n\n\nSince Bambi works with long-format data frames, we need to reshape our data before creating the data frame.\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 1),\n \"y\": np.tile(yy.flatten(), 1), \n \"outcome\": f.flatten()\n }\n)\n\nNow, let’s construct the model. The only notable distinction from the one-dimensional case is that we provide two unnamed arguments to the hsgp() function, representing the predictors on each dimension.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, c=1.5, m=10)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\"outcome ~ 0 + hsgp(x, y, c=1.5, m=10)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x, y, c=1.5, m=10)\": \"hsgp\"})\nmodel\n\n Formula: outcome ~ 0 + hsgp(x, y, c=1.5, m=10)\n Family: gaussian\n Link: mu = identity\n Observations: 144\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, c=1.5, m=10)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\nThe parameters c and m of the HSGP aproximation are specific to each dimension, and can have different values for each. However, as we are passing scalars instead of sequences, Bambi will internally recycle them, causing the HSGP approximation to use the same values of c and m for both dimensions.\nLet’s build the internal PyMC model and create a graph to have a visual representation of the relationships between the model parameters.\n\nmodel.build()\nmodel.graph()\n\n\n\n\nAnd finally, we quickly fit the model and show a traceplot to explore the posterior and spot any issues with the sampler.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\n/home/tomas/anaconda3/envs/bambi_hsgp/lib/python3.10/site-packages/pymc/sampling/jax.py:39: UserWarning: This module is experimental.\n warnings.warn(\"This module is experimental.\")\n\n\nCompiling...\n\n\nNo GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n\n\nCompilation time = 0:00:02.522713\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:23.313351\nTransforming variables...\nTransformation time = 0:00:00.628279\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nWe don’t see any divergences. However, the autocorrelation in the chains for the covariance function parameters, along with the insufficient mixing, indicates that there may be an issue with the prior specification of the model.\nSince the goal of the notebook is to simply show what features Bambi supports and how to use them, we won’t further investigate these issues. However, such posteriors shouldn’t be considered in any serious application.\nFrom now on, the notebook will follow the same structure as the one already shown, which consists of\n\nData simulation with some specific settings\nCreation of the Bambi model\nBuilding of the internal PyMC model and visualization of the graph\nModel fit and inspection of the traceplot\n\n\n\n\nIn this scenario we have multiple groups that share the same covariance function.\n\nrng = np.random.default_rng(123)\n\nell = 2\ncov = 1.2 * pm.gp.cov.ExpQuad(2, ls=ell)\nK = cov(X).eval()\nmu = np.zeros(X.shape[0])\n\nf = rng.multivariate_normal(mu, K, 3)\n\nfig, axes = plt.subplots(1, 3, figsize=(12, 4))\nfor i, ax in enumerate(axes):\n ax.scatter(xx, yy, c=f[i], s=320, marker=\"s\")\n ax.grid(False)\n ax.set_title(f\"Group {i}\")\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 3),\n \"y\": np.tile(yy.flatten(), 3),\n \"group\": np.repeat(list(\"ABC\"), 12 * 12),\n \"outcome\": f.flatten()\n }\n)\n\nNotice we don’t modify anything substantial in the call to hsgp() for now.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, by=group, c=1.5, m=10)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\"outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x, y, by=group, c=1.5, m=10)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10)\n Family: gaussian\n Link: mu = identity\n Observations: 432\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, by=group, c=1.5, m=10)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:02.721842\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:02:17.782596\nTransforming variables...\nTransformation time = 0:00:00.838094\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nWhile we have three groups, we only have one hsgp_sigma and one hsgp_ell for all groups. This is because, by default, the HSGP contributions by groups use the same instance of the covariance function.\n\n\n\nAgain we have multiple groups. But this time, each group has specific values for the amplitude and the lengthscale.\n\nrng = np.random.default_rng(12)\n\nsigmas = [1.2, 1.5, 1.8]\nells = [1.5, 2, 3]\n\nsamples = []\nfor sigma, ell in zip(sigmas, ells):\n cov = sigma * pm.gp.cov.ExpQuad(2, ls=ell)\n K = cov(X).eval()\n mu = np.zeros(X.shape[0])\n samples.append(rng.multivariate_normal(mu, K))\n\nf = np.stack(samples)\nfig, axes = plt.subplots(1, 3, figsize=(12, 4))\nfor i, ax in enumerate(axes):\n ax.scatter(xx, yy, c=f[i], s=320, marker=\"s\")\n ax.grid(False)\n ax.set_title(f\"Group {i}\")\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 3),\n \"y\": np.tile(yy.flatten(), 3),\n \"group\": np.repeat(list(\"ABC\"), 12 * 12),\n \"outcome\": f.flatten()\n }\n)\n\nIn situations like this, we can tell Bambi not to use the same covariance function for all the groups with share_cov=False and Bambi will create a separate instance for each group, resulting in group specific estimates of the amplitude and the lengthscale.\nNotice, however, we’re still using the same kind of covariance function, which in this case is ExpQuad.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\n \"outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\", \n data, \n priors=priors\n)\nmodel.set_alias({\"hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\n Family: gaussian\n Link: mu = identity\n Observations: 432\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\nSee the all the HSGP related parameters gained the new dimension hsgp_by.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:04.491697\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:02:35.274256\nTransforming variables...\nTransformation time = 0:00:00.801181\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nUnlike the previous case, now there are three hsgp_sigma and three hsgp_ell parameters, one per group. We can see them in different colors in the visualization.\n\n\n\n\nIn this second part we repeat exactly the same that we did for the isotropic case. First, we start with a single group. Then, we continue with multiple groups that share the covariance function. And finally, multiple groups with different covariance functions. The main difference is that we use iso=False, which asks to use an anisotropic GP.\n\n\n\nrng = np.random.default_rng(1234)\n\nell = [2, 0.9]\ncov = 1.2 * pm.gp.cov.ExpQuad(2, ls=ell)\nK = cov(X).eval()\nmu = np.zeros(X.shape[0])\n\nf = rng.multivariate_normal(mu, K)\n\nfig, ax = plt.subplots(figsize = (4.5, 4.5))\nax.scatter(xx, yy, c=f, s=900, marker=\"s\");\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 1),\n \"y\": np.tile(yy.flatten(), 1), \n \"outcome\": f.flatten()\n }\n)\n\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, c=1.5, m=10, iso=False)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\"outcome ~ 0 + hsgp(x, y, c=1.5, m=10, iso=False)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x, y, c=1.5, m=10, iso=False)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, c=1.5, m=10, iso=False)\n Family: gaussian\n Link: mu = identity\n Observations: 144\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, c=1.5, m=10, iso=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\nAlthough there is only one group in this case, the graph includes a hsgp_var dimension. This dimension represents the variables in the HSGP component, indicating that there is one lengthscale parameter per variable.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:02.320646\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:06.159032\nTransforming variables...\nTransformation time = 0:00:00.173091\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\n\n\n\n\nrng = np.random.default_rng(123)\n\nell = [2, 0.9]\ncov = 1.2 * pm.gp.cov.ExpQuad(2, ls=ell)\nK = cov(X).eval()\nmu = np.zeros(X.shape[0])\n\nf = rng.multivariate_normal(mu, K, 3)\n\nfig, axes = plt.subplots(1, 3, figsize=(12, 4))\nfor i, ax in enumerate(axes):\n ax.scatter(xx, yy, c=f[i], s=320, marker=\"s\")\n ax.grid(False)\n ax.set_title(f\"Group {i}\")\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 3),\n \"y\": np.tile(yy.flatten(), 3),\n \"group\": np.repeat(list(\"ABC\"), 12 * 12),\n \"outcome\": f.flatten()\n }\n)\n\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, by=group, c=1.5, m=10, iso=False)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\"outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, iso=False)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x, y, by=group, c=1.5, m=10, iso=False)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, iso=False)\n Family: gaussian\n Link: mu = identity\n Observations: 432\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, by=group, c=1.5, m=10, iso=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:02.464203\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:17.674547\nTransforming variables...\nTransformation time = 0:00:00.249682\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\n\n\n\n\nrng = np.random.default_rng(12)\n\nsigmas = [1.2, 1.5, 1.8]\nells = [[1.5, 0.8], [2, 1.5], [3, 1]]\n\nsamples = []\nfor sigma, ell in zip(sigmas, ells):\n cov = sigma * pm.gp.cov.ExpQuad(2, ls=ell)\n K = cov(X).eval()\n mu = np.zeros(X.shape[0])\n samples.append(rng.multivariate_normal(mu, K))\n\nf = np.stack(samples)\nfig, axes = plt.subplots(1, 3, figsize=(12, 4))\nfor i, ax in enumerate(axes):\n ax.scatter(xx, yy, c=f[i], s=320, marker=\"s\")\n ax.grid(False)\n ax.set_title(f\"Group {i}\")\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 3),\n \"y\": np.tile(yy.flatten(), 3),\n \"group\": np.repeat(list(\"ABC\"), 12 * 12),\n \"outcome\": f.flatten()\n }\n)\n\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\n \"outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\", \n data, \n priors=priors\n)\nmodel.set_alias({\"hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\n Family: gaussian\n Link: mu = identity\n Observations: 432\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:03.955870\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:20.713181\nTransforming variables...\nTransformation time = 0:00:00.513813\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\n\n\n\n\nFor this final demonstration we’re going to use a simulated dataset where the outcome is a count variable. For the predictors, we have the location in terms of the latitude and longitude, as well as other variables such as the year of the measurement, the site where the measure was made, and one continuous predictor.\n\ndata = pd.read_csv(\"data/poisson_data.csv\")\ndata[\"Year\"] = pd.Categorical(data[\"Year\"])\nprint(data.shape)\ndata.head()\n\n(100, 6)\n\n\n\n\n\n\n \n \n \n Year\n Count\n Site\n Lat\n Lon\n X1\n \n \n \n \n 0\n 2015\n 4\n Site1\n 47.559880\n 7.216754\n 3.316140\n \n \n 1\n 2016\n 0\n Site1\n 47.257079\n 7.135390\n 2.249612\n \n \n 2\n 2015\n 0\n Site1\n 47.061967\n 7.804383\n 2.835283\n \n \n 3\n 2016\n 0\n Site1\n 47.385533\n 7.433145\n 2.776692\n \n \n 4\n 2015\n 1\n Site1\n 47.034987\n 7.434643\n 2.295769\n \n \n\n\n\n\nWe can visualize the outcome variable by location and year.\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 4))\nfor i, (ax, year) in enumerate(zip(axes, [2015, 2016])):\n mask = data[\"Year\"] == year\n x = data.loc[mask, \"Lat\"]\n y = data.loc[mask, \"Lon\"]\n count = data.loc[mask, \"Count\"]\n ax.scatter(x, y, c=count, s=30, marker=\"s\")\n ax.set_title(f\"Year {year}\")\n\n\n\n\nThere’s not much we can conclude from here but it’s not a problem. The most relevant part of the example is not the data itself, but how to use Bambi to include GP components in a complex model.\nIt’s very easy to create a model that uses both regular common and group-specific predictors as well as a GP contribution term. We just add them to the model formula, treat hsgp() as any other call, and that’s it!\nBelow we have common effects for the Year, the interaction between X1 and Year, and group-specific intercepts by Site. Finally, we add hsgp() as any other call.\n\nformula = \"Count ~ 0 + Year + X1:Year + (1|Site) + hsgp(Lon, Lat, by=Year, m=5, c=1.5)\"\nmodel = bmb.Model(formula, data, family=\"poisson\")\nmodel\n\n Formula: Count ~ 0 + Year + X1:Year + (1|Site) + hsgp(Lon, Lat, by=Year, m=5, c=1.5)\n Family: poisson\n Link: mu = log\n Observations: 100\n Priors: \n target = mu\n Common-level effects\n Year ~ Normal(mu: [0. 0.], sigma: [5. 5.])\n X1:Year ~ Normal(mu: [0. 0.], sigma: [1.5693 1.4766])\n \n Group-level effects\n 1|Site ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 5.3683))\n \n HSGP contributions\n hsgp(Lon, Lat, by=Year, m=5, c=1.5)\n cov: ExpQuad\n sigma ~ Exponential(lam: 1.0)\n ell ~ InverseGamma(alpha: 3.0, beta: 2.0)\n\n\nLet’s use an alias to make the graph representation more readable.\n\nmodel.set_alias({\"hsgp(Lon, Lat, by=Year, m=5, c=1.5)\": \"gp\"})\nmodel.build()\nmodel.graph()\n\n\n\n\nAnd finally, let’s fit the model.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.99)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:04.433012\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:09.698066\nTransforming variables...\nTransformation time = 0:00:00.668909\n15\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"gp_sigma\", \"gp_ell\", \"gp_weights\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nNotice the posteriors for the gp_weights are all centered at zero. This is a symptom of the absence of any spatial effect.\n\naz.plot_trace(\n idata, \n var_names=[\"Year\", \"X1:Year\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"1|Site\", \"1|Site_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);" + "text": "This article demonstrates how to use Bambi with Gaussian Processes with 2 dimensional predictors. Bambi supports Gaussian Processes through the low-rank approximation known as Hilbert Space Gaussian Processes. For references see Hilbert Space Methods for Reduced-Rank Gaussian Process Regression and Practical Hilbert Space Approximate Bayesian Gaussian Processes for Probabilistic Programming.\nIf you prefer a video format, have a look at Introduction to Hilbert Space GPs in PyMC given by Bill Engels.\n\nimport arviz as az\nimport bambi as bmb\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport pymc as pm\n\nThe goal of this notebook is to showcase Bambi’s support for Gaussian Processes on two-dimensional data using the HSGP approximation.\nTo achieve this, we begin by creating a matrix of coordinates that will serve as the locations where we measure the values of a continuous response variable.\n\nx1 = np.linspace(0, 10, 12)\nx2 = np.linspace(0, 10, 12)\nxx, yy = np.meshgrid(x1, x2)\nX = np.column_stack([xx.flatten(), yy.flatten()])\nX.shape\n\n(144, 2)\n\n\n\n\nIn modeling multi-dimensional data with a Gaussian Process, we must choose between using an isotropic or an anisotropic Gaussian Process. An isotropic GP applies the same degree of smoothing to all predictors and is rotationally invariant. On the other hand, an anisotropic GP assigns different degrees of smoothing to each predictor and is not rotationally invariant.\nFurthermore, as the hsgp() function allows for the creation of separate GP contribution terms for the levels of a categorical variable through its by argument, we also examine both single-group and multiple-group scenarios.\n\n\nWe create a covariance kernel using ExpQuad from the gp submodule in PyMC. Note that the lengthscale and amplitude for both dimensions are 2 and 1.2, respectively. Then, we simply use NumPy to get a random draw from the 144-dimensional multivariate normal distribution.\n\nrng = np.random.default_rng(1234)\n\nell = 2\ncov = 1.2 * pm.gp.cov.ExpQuad(2, ls=ell)\nK = cov(X).eval()\nmu = np.zeros(X.shape[0])\nprint(mu.shape, K.shape)\n\nf = rng.multivariate_normal(mu, K)\n\nfig, ax = plt.subplots()\nax.scatter(xx, yy, c=f, s=900, marker=\"s\");\n\n(144,) (144, 144)\n\n\n\n\n\nSince Bambi works with long-format data frames, we need to reshape our data before creating the data frame.\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 1),\n \"y\": np.tile(yy.flatten(), 1), \n \"outcome\": f.flatten()\n }\n)\n\nNow, let’s construct the model. The only notable distinction from the one-dimensional case is that we provide two unnamed arguments to the hsgp() function, representing the predictors on each dimension.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, c=1.5, m=10)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\"outcome ~ 0 + hsgp(x, y, c=1.5, m=10)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x, y, c=1.5, m=10)\": \"hsgp\"})\nmodel\n\n Formula: outcome ~ 0 + hsgp(x, y, c=1.5, m=10)\n Family: gaussian\n Link: mu = identity\n Observations: 144\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, c=1.5, m=10)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\nThe parameters c and m of the HSGP aproximation are specific to each dimension, and can have different values for each. However, as we are passing scalars instead of sequences, Bambi will internally recycle them, causing the HSGP approximation to use the same values of c and m for both dimensions.\nLet’s build the internal PyMC model and create a graph to have a visual representation of the relationships between the model parameters.\n\nmodel.build()\nmodel.graph()\n\n\n\n\nAnd finally, we quickly fit the model and show a traceplot to explore the posterior and spot any issues with the sampler.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\n/home/tomas/anaconda3/envs/bambi_hsgp/lib/python3.10/site-packages/pymc/sampling/jax.py:39: UserWarning: This module is experimental.\n warnings.warn(\"This module is experimental.\")\n\n\nCompiling...\n\n\nNo GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n\n\nCompilation time = 0:00:02.522713\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:23.313351\nTransforming variables...\nTransformation time = 0:00:00.628279\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nWe don’t see any divergences. However, the autocorrelation in the chains for the covariance function parameters, along with the insufficient mixing, indicates that there may be an issue with the prior specification of the model.\nSince the goal of the notebook is to simply show what features Bambi supports and how to use them, we won’t further investigate these issues. However, such posteriors shouldn’t be considered in any serious application.\nFrom now on, the notebook will follow the same structure as the one already shown, which consists of\n\nData simulation with some specific settings\nCreation of the Bambi model\nBuilding of the internal PyMC model and visualization of the graph\nModel fit and inspection of the traceplot\n\n\n\n\nIn this scenario we have multiple groups that share the same covariance function.\n\nrng = np.random.default_rng(123)\n\nell = 2\ncov = 1.2 * pm.gp.cov.ExpQuad(2, ls=ell)\nK = cov(X).eval()\nmu = np.zeros(X.shape[0])\n\nf = rng.multivariate_normal(mu, K, 3)\n\nfig, axes = plt.subplots(1, 3, figsize=(12, 4))\nfor i, ax in enumerate(axes):\n ax.scatter(xx, yy, c=f[i], s=320, marker=\"s\")\n ax.grid(False)\n ax.set_title(f\"Group {i}\")\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 3),\n \"y\": np.tile(yy.flatten(), 3),\n \"group\": np.repeat(list(\"ABC\"), 12 * 12),\n \"outcome\": f.flatten()\n }\n)\n\nNotice we don’t modify anything substantial in the call to hsgp() for now.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, by=group, c=1.5, m=10)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\"outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x, y, by=group, c=1.5, m=10)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10)\n Family: gaussian\n Link: mu = identity\n Observations: 432\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, by=group, c=1.5, m=10)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:02.721842\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:02:17.782596\nTransforming variables...\nTransformation time = 0:00:00.838094\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nWhile we have three groups, we only have one hsgp_sigma and one hsgp_ell for all groups. This is because, by default, the HSGP contributions by groups use the same instance of the covariance function.\n\n\n\nAgain we have multiple groups. But this time, each group has specific values for the amplitude and the lengthscale.\n\nrng = np.random.default_rng(12)\n\nsigmas = [1.2, 1.5, 1.8]\nells = [1.5, 2, 3]\n\nsamples = []\nfor sigma, ell in zip(sigmas, ells):\n cov = sigma * pm.gp.cov.ExpQuad(2, ls=ell)\n K = cov(X).eval()\n mu = np.zeros(X.shape[0])\n samples.append(rng.multivariate_normal(mu, K))\n\nf = np.stack(samples)\nfig, axes = plt.subplots(1, 3, figsize=(12, 4))\nfor i, ax in enumerate(axes):\n ax.scatter(xx, yy, c=f[i], s=320, marker=\"s\")\n ax.grid(False)\n ax.set_title(f\"Group {i}\")\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 3),\n \"y\": np.tile(yy.flatten(), 3),\n \"group\": np.repeat(list(\"ABC\"), 12 * 12),\n \"outcome\": f.flatten()\n }\n)\n\nIn situations like this, we can tell Bambi not to use the same covariance function for all the groups with share_cov=False and Bambi will create a separate instance for each group, resulting in group specific estimates of the amplitude and the lengthscale.\nNotice, however, we’re still using the same kind of covariance function, which in this case is ExpQuad.\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\n \"outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\", \n data, \n priors=priors\n)\nmodel.set_alias({\"hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\n Family: gaussian\n Link: mu = identity\n Observations: 432\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, by=group, c=1.5, m=10, share_cov=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\nSee the all the HSGP related parameters gained the new dimension hsgp_by.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:04.491697\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:02:35.274256\nTransforming variables...\nTransformation time = 0:00:00.801181\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nUnlike the previous case, now there are three hsgp_sigma and three hsgp_ell parameters, one per group. We can see them in different colors in the visualization.\n\n\n\n\nIn this second part we repeat exactly the same that we did for the isotropic case. First, we start with a single group. Then, we continue with multiple groups that share the covariance function. And finally, multiple groups with different covariance functions. The main difference is that we use iso=False, which asks to use an anisotropic GP.\n\n\n\nrng = np.random.default_rng(1234)\n\nell = [2, 0.9]\ncov = 1.2 * pm.gp.cov.ExpQuad(2, ls=ell)\nK = cov(X).eval()\nmu = np.zeros(X.shape[0])\n\nf = rng.multivariate_normal(mu, K)\n\nfig, ax = plt.subplots(figsize = (4.5, 4.5))\nax.scatter(xx, yy, c=f, s=900, marker=\"s\");\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 1),\n \"y\": np.tile(yy.flatten(), 1), \n \"outcome\": f.flatten()\n }\n)\n\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, c=1.5, m=10, iso=False)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\"outcome ~ 0 + hsgp(x, y, c=1.5, m=10, iso=False)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x, y, c=1.5, m=10, iso=False)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, c=1.5, m=10, iso=False)\n Family: gaussian\n Link: mu = identity\n Observations: 144\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, c=1.5, m=10, iso=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\nAlthough there is only one group in this case, the graph includes a hsgp_var dimension. This dimension represents the variables in the HSGP component, indicating that there is one lengthscale parameter per variable.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:02.320646\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:06.159032\nTransforming variables...\nTransformation time = 0:00:00.173091\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\n\n\n\n\nrng = np.random.default_rng(123)\n\nell = [2, 0.9]\ncov = 1.2 * pm.gp.cov.ExpQuad(2, ls=ell)\nK = cov(X).eval()\nmu = np.zeros(X.shape[0])\n\nf = rng.multivariate_normal(mu, K, 3)\n\nfig, axes = plt.subplots(1, 3, figsize=(12, 4))\nfor i, ax in enumerate(axes):\n ax.scatter(xx, yy, c=f[i], s=320, marker=\"s\")\n ax.grid(False)\n ax.set_title(f\"Group {i}\")\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 3),\n \"y\": np.tile(yy.flatten(), 3),\n \"group\": np.repeat(list(\"ABC\"), 12 * 12),\n \"outcome\": f.flatten()\n }\n)\n\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, by=group, c=1.5, m=10, iso=False)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\"outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, iso=False)\", data, priors=priors)\nmodel.set_alias({\"hsgp(x, y, by=group, c=1.5, m=10, iso=False)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, iso=False)\n Family: gaussian\n Link: mu = identity\n Observations: 432\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, by=group, c=1.5, m=10, iso=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:02.464203\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:17.674547\nTransforming variables...\nTransformation time = 0:00:00.249682\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\n\n\n\n\nrng = np.random.default_rng(12)\n\nsigmas = [1.2, 1.5, 1.8]\nells = [[1.5, 0.8], [2, 1.5], [3, 1]]\n\nsamples = []\nfor sigma, ell in zip(sigmas, ells):\n cov = sigma * pm.gp.cov.ExpQuad(2, ls=ell)\n K = cov(X).eval()\n mu = np.zeros(X.shape[0])\n samples.append(rng.multivariate_normal(mu, K))\n\nf = np.stack(samples)\nfig, axes = plt.subplots(1, 3, figsize=(12, 4))\nfor i, ax in enumerate(axes):\n ax.scatter(xx, yy, c=f[i], s=320, marker=\"s\")\n ax.grid(False)\n ax.set_title(f\"Group {i}\")\n\n\n\n\n\ndata = pd.DataFrame(\n {\n \"x\": np.tile(xx.flatten(), 3),\n \"y\": np.tile(yy.flatten(), 3),\n \"group\": np.repeat(list(\"ABC\"), 12 * 12),\n \"outcome\": f.flatten()\n }\n)\n\n\nprior_hsgp = {\n \"sigma\": bmb.Prior(\"Exponential\", lam=3),\n \"ell\": bmb.Prior(\"InverseGamma\", mu=2, sigma=0.2),\n}\npriors = {\n \"hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\": prior_hsgp, \n \"sigma\": bmb.Prior(\"HalfNormal\", sigma=2)\n}\nmodel = bmb.Model(\n \"outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\", \n data, \n priors=priors\n)\nmodel.set_alias({\"hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\": \"hsgp\"})\nprint(model)\nmodel.build()\nmodel.graph()\n\n Formula: outcome ~ 0 + hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\n Family: gaussian\n Link: mu = identity\n Observations: 432\n Priors: \n target = mu\n HSGP contributions\n hsgp(x, y, by=group, c=1.5, m=10, iso=False, share_cov=False)\n cov: ExpQuad\n sigma ~ Exponential(lam: 3.0)\n ell ~ InverseGamma(mu: 2.0, sigma: 0.2)\n \n Auxiliary parameters\n outcome_sigma ~ HalfNormal(sigma: 2.0)\n\n\n\n\n\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.9)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:03.955870\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:20.713181\nTransforming variables...\nTransformation time = 0:00:00.513813\n0\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"hsgp_sigma\", \"hsgp_ell\", \"outcome_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\n\n\n\n\nFor this final demonstration we’re going to use a simulated dataset where the outcome is a count variable. For the predictors, we have the location in terms of the latitude and longitude, as well as other variables such as the year of the measurement, the site where the measure was made, and one continuous predictor.\n\ndata = pd.read_csv(\"data/poisson_data.csv\")\ndata[\"Year\"] = pd.Categorical(data[\"Year\"])\nprint(data.shape)\ndata.head()\n\n(100, 6)\n\n\n\n\n\n\n \n \n \n Year\n Count\n Site\n Lat\n Lon\n X1\n \n \n \n \n 0\n 2015\n 4\n Site1\n 47.559880\n 7.216754\n 3.316140\n \n \n 1\n 2016\n 0\n Site1\n 47.257079\n 7.135390\n 2.249612\n \n \n 2\n 2015\n 0\n Site1\n 47.061967\n 7.804383\n 2.835283\n \n \n 3\n 2016\n 0\n Site1\n 47.385533\n 7.433145\n 2.776692\n \n \n 4\n 2015\n 1\n Site1\n 47.034987\n 7.434643\n 2.295769\n \n \n\n\n\n\nWe can visualize the outcome variable by location and year.\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 4))\nfor i, (ax, year) in enumerate(zip(axes, [2015, 2016])):\n mask = data[\"Year\"] == year\n x = data.loc[mask, \"Lat\"]\n y = data.loc[mask, \"Lon\"]\n count = data.loc[mask, \"Count\"]\n ax.scatter(x, y, c=count, s=30, marker=\"s\")\n ax.set_title(f\"Year {year}\")\n\n\n\n\nThere’s not much we can conclude from here but it’s not a problem. The most relevant part of the example is not the data itself, but how to use Bambi to include GP components in a complex model.\nIt’s very easy to create a model that uses both regular common and group-specific predictors as well as a GP contribution term. We just add them to the model formula, treat hsgp() as any other call, and that’s it!\nBelow we have common effects for the Year, the interaction between X1 and Year, and group-specific intercepts by Site. Finally, we add hsgp() as any other call.\n\nformula = \"Count ~ 0 + Year + X1:Year + (1|Site) + hsgp(Lon, Lat, by=Year, m=5, c=1.5)\"\nmodel = bmb.Model(formula, data, family=\"poisson\")\nmodel\n\n Formula: Count ~ 0 + Year + X1:Year + (1|Site) + hsgp(Lon, Lat, by=Year, m=5, c=1.5)\n Family: poisson\n Link: mu = log\n Observations: 100\n Priors: \n target = mu\n Common-level effects\n Year ~ Normal(mu: [0. 0.], sigma: [5. 5.])\n X1:Year ~ Normal(mu: [0. 0.], sigma: [1.5693 1.4766])\n \n Group-level effects\n 1|Site ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 5.3683))\n \n HSGP contributions\n hsgp(Lon, Lat, by=Year, m=5, c=1.5)\n cov: ExpQuad\n sigma ~ Exponential(lam: 1.0)\n ell ~ InverseGamma(alpha: 3.0, beta: 2.0)\n\n\nLet’s use an alias to make the graph representation more readable.\n\nmodel.set_alias({\"hsgp(Lon, Lat, by=Year, m=5, c=1.5)\": \"gp\"})\nmodel.build()\nmodel.graph()\n\n\n\n\nAnd finally, let’s fit the model.\n\nidata = model.fit(inference_method=\"nuts_numpyro\", target_accept=0.99)\nprint(idata.sample_stats.diverging.sum().item())\n\nCompiling...\nCompilation time = 0:00:04.433012\nSampling...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSampling time = 0:00:09.698066\nTransforming variables...\nTransformation time = 0:00:00.668909\n15\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"gp_sigma\", \"gp_ell\", \"gp_weights\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\nNotice the posteriors for the gp_weights are all centered at zero. This is a symptom of the absence of any spatial effect.\n\naz.plot_trace(\n idata, \n var_names=[\"Year\", \"X1:Year\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);\n\n\n\n\n\naz.plot_trace(\n idata, \n var_names=[\"1|Site\", \"1|Site_sigma\"], \n backend_kwargs={\"layout\": \"constrained\"}\n);" }, { "objectID": "notebooks/logistic_regression.html",