Interpret Advanced Usage

The interpret module provides a set of helper functions to aid the user in more advanced and complex analysis not covered within the comparisons, predictions, and slopes functions. These helper functions are data_grid and select_draws. The data_grid can be used to create a pairwise grid of data points for the user to pass to model.predict. Subsequently, select_draws is used to select the draws from the posterior (or posterior predictive) group of the InferenceData object returned by the predict method that correspond to the data points that “produced” that draw.

With access to the appropriately indexed draws, and data used to generate those draws, it enables for more complex analysis such as cross-comparisons and the choice of which model parameter to compute a quantity of interest for; among others. Additionally, the user has more control over the data passed to model.predict. Below, it will be demonstrated how to use these helper functions. First, to reproduce the results from the standard interpret API, and then to compute cross-comparisons.

import warnings

import arviz as az
import numpy as np
import pandas as pd

import bambi as bmb

from bambi.interpret.helpers import data_grid, select_draws

warnings.simplefilter(action='ignore', category=FutureWarning)

Zero Inflated Poisson

We will adopt the zero inflated Poisson (ZIP) model from the comparisons documentation to demonstrate the helper functions introduced above.

The ZIP model will be used to predict how many fish are caught by visitors at a state park using survey data. Many visitors catch zero fish, either because they did not fish at all, or because they were unlucky. We would like to explicitly model this bimodal behavior (zero versus non-zero) using a Zero Inflated Poisson model, and to compare how different inputs of interest \(w\) and other covariate values \(c\) are associated with the number of fish caught. The dataset contains data on 250 groups that went to a state park to fish. Each group was questioned about how many fish they caught (count), how many children were in the group (child), how many people were in the group (persons), if they used a live bait and whether or not they brought a camper to the park (camper).

fish_data = pd.read_csv("https://stats.idre.ucla.edu/stat/data/fish.csv")
cols = ["count", "livebait", "camper", "persons", "child"]
fish_data = fish_data[cols]
fish_data["child"] = fish_data["child"].astype(np.int8)
fish_data["persons"] = fish_data["persons"].astype(np.int8)
fish_data["livebait"] = pd.Categorical(fish_data["livebait"])
fish_data["camper"] = pd.Categorical(fish_data["camper"])

fish_model = bmb.Model(
    "count ~ livebait + camper + persons + child", 
    fish_data, 
    family='zero_inflated_poisson'
)

fish_idata = fish_model.fit(random_seed=1234)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [psi, Intercept, livebait, camper, persons, child]

Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 8 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics

Create a grid of data

data_grid allows you to create a pairwise grid, also known as a cross-join or cartesian product, of data using the covariates passed to the conditional and the optional variable parameter. Covariates not passed to conditional, but are terms in the Bambi model, are set to typical values (e.g., mean or mode). If you are coming from R, this function is partially inspired from the data_grid function in {modelr}.

There are two ways to create a pairwise grid of data:

1. user-provided values are passed as a dictionary to conditional where the keys are the names of the covariates and the values are the values to use in the grid.
2. a list of covariates where the elements are the names of the covariates to use in the grid. As only the names of the covariates were passed, default values are computed to construct the grid.

Any unspecified covariates, i.e., covariates not passed to conditional but are terms in the Bambi model, are set to their “typical” values such as mean or mode depending on the data type of the covariate.

User-provided values

To construct a pairwise grid of data for specific covariate values, pass a dictionary to conditional. The values of the dictionary can be of type int, float, list, or np.ndarray.

conditional = {
    "camper": np.array([0, 1]),
    "persons": np.arange(1, 5, 1),
    "child": np.array([1, 2, 3]),
}
user_passed_grid = data_grid(fish_model, conditional)
user_passed_grid.query("camper == 0")

Default computed for unspecified variable: livebait

	camper	persons	child	livebait
0	0	1	1	1
1	0	1	2	1
2	0	1	3	1
3	0	2	1	1
4	0	2	2	1
5	0	2	3	1
6	0	3	1	1
7	0	3	2	1
8	0	3	3	1
9	0	4	1	1
10	0	4	2	1
11	0	4	3	1

Subsetting by camper = 0, it can be seen that a combination of all possible pairs of values from the dictionary (including the unspecified variable livebait) results in a dataframe containing every possible combination of values from the original sets. livebait has been set to 1 as this is the mode of the unspecified categorical variable.

Default values

Alternatively, a list of covariates can be passed to conditional where the elements are the names of the covariates to use in the grid. By doing this, you are telling interpret to compute default values for these covariates. The psuedocode below outlines the logic and functions used to compute these default values:

if is_numeric_dtype(x) or is_float_dtype(x):
    values = np.linspace(np.min(x), np.max(x), 50)

elif is_integer_dtype(x):
    values = np.quantile(x, np.linspace(0, 1, 5))

elif is_categorical_dtype(x) or is_string_dtype(x) or is_object_dtype(x):
    values = np.unique(x)

conditional = ["camper", "persons", "child"]
default_grid = data_grid(fish_model, conditional)

default_grid.shape, user_passed_grid.shape

Default computed for conditional variable: camper, persons, child
Default computed for unspecified variable: livebait

((32, 4), (24, 4))

Notice how the resulting length is different between the user passed and default grid. This is due to the fact that values for child range from 0 to 3 for the default grid.

default_grid.query("camper == 0")

	camper	persons	child	livebait
0	0	1	0	1
1	0	1	1	1
2	0	1	2	1
3	0	1	3	1
4	0	2	0	1
5	0	2	1	1
6	0	2	2	1
7	0	2	3	1
8	0	3	0	1
9	0	3	1	1
10	0	3	2	1
11	0	3	3	1
12	0	4	0	1
13	0	4	1	1
14	0	4	2	1
15	0	4	3	1

Compute comparisons

To use data_grid to help generate data in computing comparisons or slopes, additional data is passed to the optional variable parameter. The name variable is an abstraction for the comparisons parameter contrast and slopes parameter wrt. If you have used any of the interpret functions, these parameter names should be familiar and the use of data_grid should be analogous to comparisons, predictions, and slopes.

variable can also be passed user-provided data (as a dictionary), or a string indicating the name of the covariate of interest. If the latter, a default value will be computed. Additionally, if an argument is passed for variable, then the effect_type needs to be passed. This is because for comparisons and slopes an epsilon value eps needs to be determined to compute the centered and finite difference, respectively. You can also pass a value for eps as a kwarg.

conditional = {
    "camper": np.array([0, 1]),
    "persons": np.arange(1, 5, 1),
    "child": np.array([1, 2, 3, 4])
}
variable = "livebait"

grid = data_grid(fish_model, conditional, variable, effect_type="comparisons")

Default computed for contrast variable: livebait

idata_grid = fish_model.predict(fish_idata, data=grid, inplace=False)

Select draws conditional on data

The second helper function to aid in more advanced analysis is select_draws. This is a function that selects the posterior or posterior predictive draws from the ArviZ InferenceData object returned by model.predict given a conditional dictionary. The conditional dictionary represents the values that correspond to that draw.

For example, if we wanted to select posterior draws where livebait = [0, 1], then all we need to do is pass a dictionary where the key is the name of the covariate and the value is the value that we want to condition on (or select). The resulting InferenceData object will contain the draws that correspond to the data points where livebait = [0, 1]. Additionally, you must pass the InferenceData object returned by model.predict, the data used to generate the predictions, and the name of the data variable data_var you would like to select from the InferenceData posterior group. If you specified to return the posterior predictive samples by passing model.predict(..., kind="pps"), you can use this group instead of the posterior group by passing group="posterior_predictive".

Below, it is demonstrated how to compute comparisons for count_mean for the contrast livebait = [0, 1] using the posterior draws.

idata_grid

arviz.InferenceData

• posterior

  <xarray.Dataset> Size: 1MB
  Dimensions:       (chain: 2, draw: 1000, camper_dim: 1, livebait_dim: 1,
                     __obs__: 64)
  Coordinates:
    * chain         (chain) int64 16B 0 1
    * draw          (draw) int64 8kB 0 1 2 3 4 5 6 ... 993 994 995 996 997 998 999
    * camper_dim    (camper_dim) <U1 4B '1'
    * livebait_dim  (livebait_dim) <U1 4B '1'
    * __obs__       (__obs__) int64 512B 0 1 2 3 4 5 6 7 ... 57 58 59 60 61 62 63
  Data variables:
      Intercept     (chain, draw) float64 16kB -2.319 -2.441 ... -2.705 -2.571
      camper        (chain, draw, camper_dim) float64 16kB 0.7341 0.6539 ... 0.763
      child         (chain, draw) float64 16kB -1.232 -1.251 ... -1.437 -1.344
      livebait      (chain, draw, livebait_dim) float64 16kB 1.509 1.634 ... 1.7
      persons       (chain, draw) float64 16kB 0.8448 0.8763 ... 0.903 0.8508
      psi           (chain, draw) float64 16kB 0.6017 0.6008 ... 0.6933 0.5425
      mu            (chain, draw, __obs__) float64 1MB 0.06677 0.3021 ... 0.1248
  Attributes:
      created_at:                  2024-05-26T21:58:43.655562+00:00
      arviz_version:               0.18.0
      inference_library:           pymc
      inference_library_version:   5.15.0+23.g19be124e
      sampling_time:               8.18968677520752
      tuning_steps:                1000
      modeling_interface:          bambi
      modeling_interface_version:  0.13.1.dev39+gb7d6a6cb

  xarray.Dataset

  • Dimensions:
    • chain: 2
    • draw: 1000
    • camper_dim: 1
    • livebait_dim: 1
    • __obs__: 64
  • Coordinates: (5)
    • chain
      (chain)
      int64
      0 1

      array([0, 1])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • camper_dim
      (camper_dim)
      <U1
      '1'

      array(['1'], dtype='<U1')

    • livebait_dim
      (livebait_dim)
      <U1
      '1'

      array(['1'], dtype='<U1')

    • __obs__
      (__obs__)
      int64
      0 1 2 3 4 5 6 ... 58 59 60 61 62 63

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59, 60, 61, 62, 63])

  • Data variables: (7)
    • Intercept
      (chain, draw)
      float64
      -2.319 -2.441 ... -2.705 -2.571

      array([[-2.31938043, -2.44096796, -2.50242028, ..., -2.41008609,
              -2.15482426, -2.20237051],
             [-1.86640536, -2.67868737, -2.6339482 , ..., -2.12677232,
              -2.70501072, -2.57088187]])

    • camper
      (chain, draw, camper_dim)
      float64
      0.7341 0.6539 ... 0.5571 0.763

      array([[[0.73412642],
              [0.65388967],
              [0.7163802 ],
              ...,
              [0.7833789 ],
              [0.65108041],
              [0.66538528]],
      
             [[0.77244764],
              [0.56742013],
              [0.5377857 ],
              ...,
              [0.76498592],
              [0.55708414],
              [0.76296319]]])

    • child
      (chain, draw)
      float64
      -1.232 -1.251 ... -1.437 -1.344

      array([[-1.2318687 , -1.25127852, -1.31269468, ..., -1.3926868 ,
              -1.25774522, -1.33447648],
             [-1.28422086, -1.34797375, -1.38221129, ..., -1.40623727,
              -1.43746916, -1.34408764]])

    • livebait
      (chain, draw, livebait_dim)
      float64
      1.509 1.634 1.678 ... 1.859 1.7

      array([[[1.50949402],
              [1.63437824],
              [1.67808319],
              ...,
              [1.63405392],
              [1.42685866],
              [1.47076256]],
      
             [[1.10221442],
              [2.23616858],
              [2.14968211],
              ...,
              [1.39522862],
              [1.85878726],
              [1.6999231 ]]])

    • persons
      (chain, draw)
      float64
      0.8448 0.8763 ... 0.903 0.8508

      array([[0.84478698, 0.87628243, 0.85830801, ..., 0.83912448, 0.84743311,
              0.84884691],
             [0.82922748, 0.78660974, 0.80622882, ..., 0.81224848, 0.90296207,
              0.85075257]])

    • psi
      (chain, draw)
      float64
      0.6017 0.6008 ... 0.6933 0.5425

      array([[0.60167917, 0.60079901, 0.56511084, ..., 0.57885181, 0.56435481,
              0.49627898],
             [0.57723187, 0.58138806, 0.57220952, ..., 0.56246952, 0.69328445,
              0.54246752]])

    • mu
      (chain, draw, __obs__)
      float64
      0.06677 0.3021 ... 0.02279 0.1248

      array([[[0.06677262, 0.30210878, 0.0194807 , ..., 0.67555671,
               0.04356153, 0.19709156],
              [0.059847  , 0.30679184, 0.01712454, ..., 0.66937831,
               0.03736344, 0.19153506],
              [0.05198464, 0.27839237, 0.01398875, ..., 0.5418268 ,
               0.02722589, 0.14580232],
              ...,
              [0.0516302 , 0.26458453, 0.0128253 , ..., 0.4429866 ,
               0.02147305, 0.11004097],
              [0.07690869, 0.32037032, 0.0218647 , ..., 0.63103849,
               0.04306726, 0.17940066],
              [0.06801683, 0.29604686, 0.01790856, ..., 0.50953105,
               0.03082271, 0.13415748]],
      
             [[0.09813622, 0.29547108, 0.0271706 , ..., 0.5900807 ,
               0.05426199, 0.16337341],
              [0.03916188, 0.36645381, 0.01017293, ..., 0.4618196 ,
               0.01282033, 0.11996507],
              [0.04035941, 0.3463697 , 0.01013113, ..., 0.41970676,
               0.01227621, 0.10535599],
              ...,
              [0.06582464, 0.26566146, 0.01613123, ..., 0.39208537,
               0.02380781, 0.09608587],
              [0.03918278, 0.25139483, 0.00930701, ..., 0.37168873,
               0.01376047, 0.08828654],
              [0.04669039, 0.25556108, 0.0121758 , ..., 0.47842863,
               0.02279398, 0.12476343]]])

  • Indexes: (5)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • camper_dim
      PandasIndex

      PandasIndex(Index(['1'], dtype='object', name='camper_dim'))

    • livebait_dim
      PandasIndex

      PandasIndex(Index(['1'], dtype='object', name='livebait_dim'))

    • __obs__
      PandasIndex

      PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59, 60, 61, 62, 63],
            dtype='int64', name='__obs__'))

  • Attributes: (8)

    created_at :

    2024-05-26T21:58:43.655562+00:00

    arviz_version :

    0.18.0

    inference_library :

    pymc

    inference_library_version :

    5.15.0+23.g19be124e

    sampling_time :

    8.18968677520752

    tuning_steps :

    1000

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.1.dev39+gb7d6a6cb

  
  
• sample_stats

  <xarray.Dataset> Size: 252kB
  Dimensions:                (chain: 2, draw: 1000)
  Coordinates:
    * chain                  (chain) int64 16B 0 1
    * draw                   (draw) int64 8kB 0 1 2 3 4 5 ... 995 996 997 998 999
  Data variables: (12/17)
      acceptance_rate        (chain, draw) float64 16kB 0.4852 0.4691 ... 0.4878
      diverging              (chain, draw) bool 2kB False False ... False False
      energy                 (chain, draw) float64 16kB 756.7 754.8 ... 755.2
      energy_error           (chain, draw) float64 16kB 0.1849 -0.1576 ... 0.6626
      index_in_trajectory    (chain, draw) int64 16kB 6 -2 -2 5 3 ... 5 -2 -5 -5 6
      largest_eigval         (chain, draw) float64 16kB nan nan nan ... nan nan
      ...                     ...
      process_time_diff      (chain, draw) float64 16kB 0.004515 ... 0.002555
      reached_max_treedepth  (chain, draw) bool 2kB False False ... False False
      smallest_eigval        (chain, draw) float64 16kB nan nan nan ... nan nan
      step_size              (chain, draw) float64 16kB 0.4493 0.4493 ... 0.2491
      step_size_bar          (chain, draw) float64 16kB 0.4092 0.4092 ... 0.3806
      tree_depth             (chain, draw) int64 16kB 4 3 3 4 3 3 ... 4 3 3 3 4 3
  Attributes:
      created_at:                  2024-05-26T21:58:43.673049+00:00
      arviz_version:               0.18.0
      inference_library:           pymc
      inference_library_version:   5.15.0+23.g19be124e
      sampling_time:               8.18968677520752
      tuning_steps:                1000
      modeling_interface:          bambi
      modeling_interface_version:  0.13.1.dev39+gb7d6a6cb

  xarray.Dataset

  • Dimensions:
    • chain: 2
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1

      array([0, 1])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

  • Data variables: (17)
    • acceptance_rate
      (chain, draw)
      float64
      0.4852 0.4691 ... 0.8185 0.4878

      array([[0.48524273, 0.46907635, 0.9900161 , ..., 0.94979748, 0.40116304,
              0.93741625],
             [0.58873199, 0.91215294, 0.49416685, ..., 0.53753997, 0.81845045,
              0.48775355]])

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • energy
      (chain, draw)
      float64
      756.7 754.8 752.2 ... 756.3 755.2

      array([[756.66437695, 754.81526077, 752.2446864 , ..., 755.57787795,
              755.37790244, 754.18463256],
             [757.69096989, 758.29295824, 756.92029865, ..., 759.27488214,
              756.33063752, 755.1645284 ]])

    • energy_error
      (chain, draw)
      float64
      0.1849 -0.1576 ... -0.5226 0.6626

      array([[ 0.18490812, -0.15761089, -0.01500881, ..., -0.61241546,
               0.04715335,  0.13670632],
             [-0.30343495,  0.1003043 , -0.00980192, ...,  0.48894455,
              -0.52259825,  0.66255443]])

    • index_in_trajectory
      (chain, draw)
      int64
      6 -2 -2 5 3 -6 ... 11 5 -2 -5 -5 6

      array([[ 6, -2, -2, ..., 13,  4, -3],
             [-4, -5, -2, ..., -5, -5,  6]])

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • lp
      (chain, draw)
      float64
      -751.7 -751.0 ... -752.3 -752.8

      array([[-751.74419455, -751.04504395, -750.54322651, ..., -750.83344139,
              -751.48798   , -753.41881499],
             [-754.58377843, -754.16645365, -753.40879282, ..., -753.23969952,
              -752.34621461, -752.81312555]])

    • max_energy_error
      (chain, draw)
      float64
      1.793 2.415 ... 0.7354 1.872

      array([[ 1.79320886,  2.41509386, -0.05842541, ..., -0.61241546,
               3.19523992,  0.13670632],
             [ 1.47885656,  0.22641457,  2.85988684, ...,  1.16100195,
               0.73536126,  1.87172503]])

    • n_steps
      (chain, draw)
      float64
      15.0 7.0 7.0 15.0 ... 7.0 15.0 7.0

      array([[15.,  7.,  7., ..., 15., 15.,  7.],
             [ 7.,  7.,  7., ...,  7., 15.,  7.]])

    • perf_counter_diff
      (chain, draw)
      float64
      0.004514 0.002346 ... 0.002554

      array([[0.00451421, 0.00234647, 0.00234609, ..., 0.00520036, 0.00660593,
              0.00227478],
             [0.00251641, 0.00256408, 0.00261551, ..., 0.00242392, 0.00476616,
              0.0025536 ]])

    • perf_counter_start
      (chain, draw)
      float64
      9.441e+03 9.441e+03 ... 9.444e+03

      array([[9440.65153013, 9440.65626809, 9440.65882007, ..., 9444.36672042,
              9444.37217617, 9444.37899446],
             [9440.65396041, 9440.65673697, 9440.65951651, ..., 9444.07317971,
              9444.07582181, 9444.08080841]])

    • process_time_diff
      (chain, draw)
      float64
      0.004515 0.00233 ... 0.002555

      array([[0.00451509, 0.00233032, 0.00234784, ..., 0.00520224, 0.0066074 ,
              0.00227581],
             [0.00251693, 0.00256571, 0.00261747, ..., 0.00242573, 0.00476705,
              0.00255493]])

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • step_size
      (chain, draw)
      float64
      0.4493 0.4493 ... 0.2491 0.2491

      array([[0.44934442, 0.44934442, 0.44934442, ..., 0.44934442, 0.44934442,
              0.44934442],
             [0.24907212, 0.24907212, 0.24907212, ..., 0.24907212, 0.24907212,
              0.24907212]])

    • step_size_bar
      (chain, draw)
      float64
      0.4092 0.4092 ... 0.3806 0.3806

      array([[0.40917257, 0.40917257, 0.40917257, ..., 0.40917257, 0.40917257,
              0.40917257],
             [0.38056482, 0.38056482, 0.38056482, ..., 0.38056482, 0.38056482,
              0.38056482]])

    • tree_depth
      (chain, draw)
      int64
      4 3 3 4 3 3 3 4 ... 4 3 4 3 3 3 4 3

      array([[4, 3, 3, ..., 4, 4, 3],
             [3, 3, 3, ..., 3, 4, 3]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

  • Attributes: (8)

    created_at :

    2024-05-26T21:58:43.673049+00:00

    arviz_version :

    0.18.0

    inference_library :

    pymc

    inference_library_version :

    5.15.0+23.g19be124e

    sampling_time :

    8.18968677520752

    tuning_steps :

    1000

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.1.dev39+gb7d6a6cb

  
  
• observed_data

  <xarray.Dataset> Size: 4kB
  Dimensions:  (__obs__: 250)
  Coordinates:
    * __obs__  (__obs__) int64 2kB 0 1 2 3 4 5 6 7 ... 243 244 245 246 247 248 249
  Data variables:
      count    (__obs__) int64 2kB 0 0 0 0 1 0 0 0 0 1 0 ... 4 1 1 0 1 0 0 0 0 0 0
  Attributes:
      created_at:                  2024-05-26T21:58:43.678881+00:00
      arviz_version:               0.18.0
      inference_library:           pymc
      inference_library_version:   5.15.0+23.g19be124e
      modeling_interface:          bambi
      modeling_interface_version:  0.13.1.dev39+gb7d6a6cb

  xarray.Dataset

  • Dimensions:
    • __obs__: 250
  • Coordinates: (1)
    • __obs__
      (__obs__)
      int64
      0 1 2 3 4 5 ... 245 246 247 248 249

      array([  0,   1,   2, ..., 247, 248, 249])

  • Data variables: (1)
    • count
      (__obs__)
      int64
      0 0 0 0 1 0 0 0 ... 0 1 0 0 0 0 0 0

      array([  0,   0,   0,   0,   1,   0,   0,   0,   0,   1,   0,   0,   1,
               2,   0,   1,   0,   0,   1,   0,   1,   5,   0,   3,  30,   0,
              13,   0,   0,   0,   0,   0,  11,   5,   0,   1,   1,   7,   0,
              14,   0,  32,   0,   1,   0,   0,   0,   1,   5,   0,   1,   0,
              22,   0,  15,   0,   0,   0,   5,   4,   2,   0,   2,  32,   0,
               0,   1,   0,   0,   0,   7,   0,   0,   0,   0,   0,   0,   0,
               0,   2,   3,   1,   5,   0,   2,   1,   0,   1, 149,   0,   1,
               0,   0,   1,   0,   0,   0,   2,   2,  29,   3,   0,   0,   5,
               0,   0,   0,   0,   0,   1,   7,   1,   0,   2,   0,   2,   0,
               0,   0,   1,   0,   0,   0,   0,   0,   3,   4,   3,   3,   8,
               2,   1,   6,   0,   0,   5,   3,  31,   0,   2,   0,   0,   0,
               0,   0,   0,   6,   9,   0,   0,   0,   0,   0,   2,  15,   1,
               2,   3,   0,  65,   5,   0,   0,   0,   0,   1,   8,   0,   0,
               0,   2,   4,   5,   9,   0,   0,   0,   0,  21,   0,   6,   0,
               0,   0,   0,  16,   0,   0,   4,   2,  10,   0,   0,   0,   2,
               1,   3,   0,   0,  21,   0,   0,   2,   0,   3,   0,  38,   0,
               0,   0,   1,   3,   0,   1,   0,   0,   0,   0,   5,   0,   0,
               2,   0,   0,   0,   1,   4,   0,   0,   2,   3,   0,   0,   0,
               0,   1,   2,   0,   6,   4,   1,   1,   0,   1,   0,   0,   0,
               0,   0,   0])

  • Indexes: (1)
    • __obs__
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             240, 241, 242, 243, 244, 245, 246, 247, 248, 249],
            dtype='int64', name='__obs__', length=250))

  • Attributes: (6)

    created_at :

    2024-05-26T21:58:43.678881+00:00

    arviz_version :

    0.18.0

    inference_library :

    pymc

    inference_library_version :

    5.15.0+23.g19be124e

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.1.dev39+gb7d6a6cb

  
  

draw_1 = select_draws(idata_grid, grid, {"livebait": 0}, "mu")

draw_1 = select_draws(idata_grid, grid, {"livebait": 0}, "mu")
draw_2 = select_draws(idata_grid, grid, {"livebait": 1}, "mu")

comparison_mean = (draw_2 - draw_1).mean(("chain", "draw"))
comparison_hdi = az.hdi(draw_2 - draw_1)

comparison_df = pd.DataFrame(
    {
        "mean": comparison_mean.values,
        "hdi_low": comparison_hdi.sel(hdi="lower")["mu"].values,
        "hdi_high": comparison_hdi.sel(hdi="higher")["mu"].values,
    }
)
comparison_df.head(10)

	mean	hdi_low	hdi_high
0	0.216038	0.150384	0.294501
1	0.054384	0.029894	0.077987
2	0.013811	0.005794	0.021496
3	0.003539	0.001142	0.006083
4	0.515839	0.375190	0.670931
5	0.129765	0.078764	0.182346
6	0.032932	0.014961	0.050478
7	0.008432	0.002665	0.014138
8	1.234183	0.926056	1.553555
9	0.310259	0.195251	0.425951

We can compare this comparison with bmb.interpret.comparisons.

summary_df = bmb.interpret.comparisons(
    fish_model,
    fish_idata,
    contrast={"livebait": [0, 1]},
    conditional=conditional
)
summary_df.head(10)

	term	estimate_type	value	camper	persons	child	estimate	lower_3.0%	upper_97.0%
0	livebait	diff	(0, 1)	0	1	1	0.216038	0.150384	0.294501
1	livebait	diff	(0, 1)	0	1	2	0.054384	0.029894	0.077987
2	livebait	diff	(0, 1)	0	1	3	0.013811	0.005794	0.021496
3	livebait	diff	(0, 1)	0	1	4	0.003539	0.001142	0.006083
4	livebait	diff	(0, 1)	0	2	1	0.515839	0.375190	0.670931
5	livebait	diff	(0, 1)	0	2	2	0.129765	0.078764	0.182346
6	livebait	diff	(0, 1)	0	2	3	0.032932	0.014961	0.050478
7	livebait	diff	(0, 1)	0	2	4	0.008432	0.002665	0.014138
8	livebait	diff	(0, 1)	0	3	1	1.234183	0.926056	1.553555
9	livebait	diff	(0, 1)	0	3	2	0.310259	0.195251	0.425951

Albeit the other information in the summary_df, the columns estimate, lower_3.0%, upper_97.0% are identical.

Cross comparisons

Computing a cross-comparison is useful for when we want to compare contrasts when two (or more) predictors change at the same time. Cross-comparisons are currently not supported in the comparisons function, but we can use select_draws to compute them. For example, imagine we are interested in computing the cross-comparison between the two rows below.

summary_df.iloc[:2]

	term	estimate_type	value	camper	persons	child	estimate	lower_3.0%	upper_97.0%
0	livebait	diff	(0, 1)	0	1	1	0.216038	0.150384	0.294501
1	livebait	diff	(0, 1)	0	1	2	0.054384	0.029894	0.077987

The cross-comparison amounts to first computing the comparison for row 0, given below, and can be verified by looking at the estimate in summary_df.

cond_10 = {
    "camper": 0,
    "persons": 1,
    "child": 1,
    "livebait": 0 
}

cond_11 = {
    "camper": 0,
    "persons": 1,
    "child": 1,
    "livebait": 1
}

draws_10 = select_draws(idata_grid, grid, cond_10, "mu")
draws_11 = select_draws(idata_grid, grid, cond_11, "mu")

(draws_11 - draws_10).mean(("chain", "draw")).item()

0.21603768204145285

Next, we need to compute the comparison for row 1.

cond_20 = {
    "camper": 0,
    "persons": 1,
    "child": 2,
    "livebait": 0
}

cond_21 = {
    "camper": 0,
    "persons": 1,
    "child": 2,
    "livebait": 1
}

draws_20 = select_draws(idata_grid, grid, cond_20, "mu")
draws_21 = select_draws(idata_grid, grid, cond_21, "mu")

(draws_21 - draws_20).mean(("chain", "draw")).item()

0.0543836038546379

Next, we compute the “first level” comparisons (diff_1 and diff_2). Subsequently, we compute the difference between these two differences to obtain the cross-comparison.

diff_1 = (draws_11 - draws_10)
diff_2 = (draws_21 - draws_20)

cross_comparison = (diff_2 - diff_1).mean(("chain", "draw")).item()
cross_comparison

-0.16165407818681496

To verify this is correct, we can check by performing the cross-comparison directly on the summary_df.

summary_df.iloc[1]["estimate"] - summary_df.iloc[0]["estimate"]

-0.16165407818681496

Summary

In this notebook, the interpret helper functions data_grid and select_draws were introduced and it was demonstrated how they can be used to compute pairwise grids of data and cross-comparisons. With these functions, it is left to the user to generate their grids of data and quantities of interest allowing for more flexibility and control over the type of data passed to model.predict and the quantities of interest computed.

%load_ext watermark
%watermark -n -u -v -iv -w

Last updated: Sun May 26 2024

Python implementation: CPython
Python version       : 3.11.9
IPython version      : 8.24.0

numpy : 1.26.4
arviz : 0.18.0
bambi : 0.13.1.dev39+gb7d6a6cb
pandas: 2.2.2

Watermark: 2.4.3