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)WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.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_stata("http://www.stata-press.com/data/r11/fish.dta")
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 (4 chains in 4 jobs)
NUTS: [count_psi, Intercept, livebait, camper, persons, child]Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 3 seconds.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:
conditional where the keys are the names of the covariates and the values are the values to use in 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.
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.0 |
| 1 | 0 | 1 | 2 | 1.0 |
| 2 | 0 | 1 | 3 | 1.0 |
| 3 | 0 | 2 | 1 | 1.0 |
| 4 | 0 | 2 | 2 | 1.0 |
| 5 | 0 | 2 | 3 | 1.0 |
| 6 | 0 | 3 | 1 | 1.0 |
| 7 | 0 | 3 | 2 | 1.0 |
| 8 | 0 | 3 | 3 | 1.0 |
| 9 | 0 | 4 | 1 | 1.0 |
| 10 | 0 | 4 | 2 | 1.0 |
| 11 | 0 | 4 | 3 | 1.0 |
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.
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.shapeDefault 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.
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: livebaitThe 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.
<xarray.Dataset>
Dimensions: (chain: 4, draw: 1000, livebait_dim: 1, camper_dim: 1,
count_obs: 64)
Coordinates:
* chain (chain) int64 0 1 2 3
* draw (draw) int64 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999
* livebait_dim (livebait_dim) <U3 '1.0'
* camper_dim (camper_dim) <U3 '1.0'
* count_obs (count_obs) int64 0 1 2 3 4 5 6 7 ... 56 57 58 59 60 61 62 63
Data variables:
Intercept (chain, draw) float64 -2.454 -2.31 -2.91 ... -2.652 -2.887
livebait (chain, draw, livebait_dim) float64 1.629 1.58 ... 1.799 1.967
camper (chain, draw, camper_dim) float64 0.7037 0.7089 ... 0.7128
persons (chain, draw) float64 0.8707 0.8369 0.9457 ... 0.8847 0.912
child (chain, draw) float64 -1.345 -1.412 -1.418 ... -1.293 -1.573
count_psi (chain, draw) float64 0.6311 0.6201 0.6342 ... 0.6768 0.5745
count_mean (chain, draw, count_obs) float64 0.05349 0.2728 ... 0.05777
Attributes:
created_at: 2023-12-05T18:56:31.591639
arviz_version: 0.16.1
inference_library: pymc
inference_library_version: 5.8.1
sampling_time: 2.6078336238861084
tuning_steps: 1000
modeling_interface: bambi
modeling_interface_version: 0.13.0.dev0array([0, 1, 2, 3])
array([ 0, 1, 2, ..., 997, 998, 999])
array(['1.0'], dtype='<U3')
array(['1.0'], dtype='<U3')
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])array([[-2.45369792, -2.31017959, -2.91023957, ..., -2.12445232,
-2.13810127, -2.40815127],
[-2.21843503, -2.81570517, -2.81570517, ..., -3.23806838,
-2.43778314, -2.40094452],
[-2.57030767, -2.71197967, -3.13096258, ..., -2.75122519,
-2.81845915, -2.77702878],
[-2.27327526, -2.76201826, -2.72978474, ..., -2.43555099,
-2.65150277, -2.88742403]])array([[[1.62921077],
[1.57971577],
[1.85814553],
...,
[1.39320044],
[1.31938529],
[2.0303655 ]],
[[1.44441014],
[2.03180494],
[2.03180494],
...,
[2.00383874],
[2.0219552 ],
[1.77554631]],
[[1.44183823],
[2.12560898],
[1.84301056],
...,
[1.65875902],
[1.69372851],
[1.78166201]],
[[1.64438961],
[1.88379206],
[1.95128032],
...,
[1.73718918],
[1.79946112],
[1.96668979]]])array([[[0.70368097],
[0.70888183],
[0.6432408 ],
...,
[0.64457622],
[0.68265053],
[0.59371655]],
[[0.73004824],
[0.63806292],
[0.63806292],
...,
[0.92181193],
[0.44431118],
[0.52257123]],
[[0.67251782],
[0.61970006],
[0.80110109],
...,
[0.71626836],
[0.7386686 ],
[0.80126267]],
[[0.68695406],
[0.66147853],
[0.7129354 ],
...,
[0.62523222],
[0.64449304],
[0.71282881]]])array([[0.87073326, 0.83691824, 0.94570525, ..., 0.86127639, 0.85807333,
0.7827932 ],
[0.85684603, 0.88279463, 0.88279463, ..., 0.93383647, 0.82809407,
0.84387704],
[0.93868593, 0.83511835, 0.9966305 , ..., 0.92864502, 0.9583216 ,
0.9046778 ],
[0.82278573, 0.89846761, 0.85979591, ..., 0.85493891, 0.88465389,
0.91197904]])array([[-1.34537519, -1.41232517, -1.41808649, ..., -1.28364874,
-1.21311483, -1.343049 ],
[-1.36022475, -1.38360737, -1.38360737, ..., -1.34137848,
-1.37005354, -1.32761054],
[-1.42361978, -1.45390876, -1.59921868, ..., -1.57170719,
-1.60156825, -1.49976267],
[-1.33972561, -1.29366616, -1.38679881, ..., -1.29405166,
-1.29343353, -1.57284431]])array([[0.63114698, 0.62012886, 0.63423439, ..., 0.57386427, 0.56121415,
0.57771985],
[0.61291154, 0.607627 , 0.607627 , ..., 0.63015121, 0.54061116,
0.54315861],
[0.60778057, 0.68015035, 0.57385421, ..., 0.75656914, 0.74951665,
0.68320354],
[0.61202421, 0.58032736, 0.57834964, ..., 0.65743098, 0.67681655,
0.57454481]])array([[[0.05348576, 0.27276925, 0.01392994, ..., 0.50966648,
0.02602794, 0.13273855],
[0.05582204, 0.27093652, 0.01359692, ..., 0.40216834,
0.02018278, 0.09795866],
[0.03395834, 0.21773528, 0.00822393, ..., 0.41466193,
0.01566191, 0.10042158],
...,
[0.07833 , 0.31549128, 0.02169934, ..., 0.61108678,
0.04203026, 0.16928611],
[0.08264981, 0.30920293, 0.0245693 , ..., 0.70955544,
0.05638135, 0.21092947],
[0.0513851 , 0.3913936 , 0.013414 , ..., 0.50558281,
0.01732754, 0.13198164]],
[[0.06575538, 0.27876013, 0.01687303, ..., 0.49794444,
0.03014001, 0.12777409],
[0.03627894, 0.27673 , 0.00909414, ..., 0.46511021,
0.01528485, 0.11659041],
[0.03627894, 0.27673 , 0.00909414, ..., 0.46511021,
0.01528485, 0.11659041],
...
[0.03356447, 0.17630703, 0.00697101, ..., 0.25240022,
0.00997967, 0.05242108],
[0.03137619, 0.17067787, 0.00632482, ..., 0.25730828,
0.00953508, 0.05186824],
[0.03431703, 0.20383353, 0.00765898, ..., 0.34140667,
0.01282825, 0.07619621]],
[[0.061408 , 0.31796133, 0.01608383, ..., 0.51172625,
0.02588528, 0.13403007],
[0.04254398, 0.27987149, 0.01166826, ..., 0.60418447,
0.02518935, 0.16570571],
[0.03851191, 0.271035 , 0.00962312, ..., 0.45530771,
0.01616574, 0.11376952],
...,
[0.05643511, 0.32062773, 0.01547212, ..., 0.5853596 ,
0.02824695, 0.16048086],
[0.04687446, 0.28342116, 0.01285894, ..., 0.57739214,
0.02619653, 0.1583944 ],
[0.02877382, 0.2056459 , 0.00596925, ..., 0.27844848,
0.00808248, 0.05776533]]])PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))
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))PandasIndex(Index(['1.0'], dtype='object', name='livebait_dim'))
PandasIndex(Index(['1.0'], dtype='object', name='camper_dim'))
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='count_obs'))<xarray.Dataset>
Dimensions: (chain: 4, draw: 1000)
Coordinates:
* chain (chain) int64 0 1 2 3
* draw (draw) int64 0 1 2 3 4 5 ... 994 995 996 997 998 999
Data variables: (12/17)
max_energy_error (chain, draw) float64 0.595 1.645 ... -0.2578 0.4279
reached_max_treedepth (chain, draw) bool False False False ... False False
process_time_diff (chain, draw) float64 0.000885 0.000874 ... 0.001667
diverging (chain, draw) bool False False False ... False False
n_steps (chain, draw) float64 7.0 7.0 7.0 ... 7.0 11.0 15.0
lp (chain, draw) float64 -750.1 -750.2 ... -751.9 -754.0
... ...
largest_eigval (chain, draw) float64 nan nan nan nan ... nan nan nan
step_size (chain, draw) float64 0.415 0.415 ... 0.3918 0.3918
smallest_eigval (chain, draw) float64 nan nan nan nan ... nan nan nan
perf_counter_start (chain, draw) float64 5.346e+04 ... 5.346e+04
acceptance_rate (chain, draw) float64 0.7512 0.5492 ... 0.9919 0.8264
index_in_trajectory (chain, draw) int64 6 4 -6 -5 -2 4 ... 7 1 4 4 -1 9
Attributes:
created_at: 2023-12-05T18:56:31.597938
arviz_version: 0.16.1
inference_library: pymc
inference_library_version: 5.8.1
sampling_time: 2.6078336238861084
tuning_steps: 1000
modeling_interface: bambi
modeling_interface_version: 0.13.0.dev0array([0, 1, 2, 3])
array([ 0, 1, 2, ..., 997, 998, 999])
array([[ 0.595031 , 1.64520146, 0.21527509, ..., 0.78075972,
5.80798519, 0.43782396],
[-0.37301973, -0.20402685, 2.5627162 , ..., 0.36741839,
3.6828103 , -0.50001872],
[ 1.27621281, -0.68221199, 0.44120051, ..., 2.43620146,
-0.08386685, -0.22408853],
[ 0.14004462, 0.99694174, -0.07864808, ..., 0.22969662,
-0.25775022, 0.42794799]])array([[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False]])array([[0.000885, 0.000874, 0.000876, ..., 0.000881, 0.000447, 0.000879],
[0.00087 , 0.000882, 0.000882, ..., 0.000867, 0.001703, 0.000869],
[0.000891, 0.000885, 0.000883, ..., 0.000876, 0.000443, 0.000866],
[0.000958, 0.000868, 0.001685, ..., 0.000881, 0.001296, 0.001667]])array([[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False]])array([[ 7., 7., 7., ..., 7., 3., 7.],
[ 7., 7., 7., ..., 7., 15., 7.],
[ 7., 7., 7., ..., 7., 3., 7.],
[ 7., 7., 15., ..., 7., 11., 15.]])array([[-750.0635191 , -750.19379111, -751.09897583, ..., -751.87190426,
-752.95441117, -754.23172742],
[-751.17173283, -750.81535676, -750.81535676, ..., -755.33966632,
-755.44844073, -752.33208106],
[-754.06157469, -753.90911474, -759.1519945 , ..., -755.85090606,
-755.67682047, -752.00709578],
[-750.85175653, -751.04878911, -750.6370331 , ..., -751.95509665,
-751.89343509, -753.95045712]])array([[ 4.37855069e-02, 2.70547404e-02, -5.88277115e-04, ...,
8.64705248e-02, 3.20765941e-01, 4.37823957e-01],
[-9.78110746e-02, -4.02245905e-02, 0.00000000e+00, ...,
2.05920633e-01, 2.91692817e-01, -2.47079734e-01],
[ 1.14647984e+00, -6.82211987e-01, 4.41200505e-01, ...,
-3.52379476e-01, 1.56742695e-03, -1.73449488e-01],
[ 1.45574755e-02, -4.02199593e-03, -5.80109339e-02, ...,
-9.46547980e-03, -1.22382092e-01, 4.15048801e-01]])array([[0.45015324, 0.45015324, 0.45015324, ..., 0.45015324, 0.45015324,
0.45015324],
[0.42774415, 0.42774415, 0.42774415, ..., 0.42774415, 0.42774415,
0.42774415],
[0.44696235, 0.44696235, 0.44696235, ..., 0.44696235, 0.44696235,
0.44696235],
[0.45405445, 0.45405445, 0.45405445, ..., 0.45405445, 0.45405445,
0.45405445]])array([[751.94493597, 752.16775838, 751.86589668, ..., 757.02411223,
757.59305462, 755.76371432],
[754.28875075, 752.65123698, 755.84607135, ..., 759.92489206,
761.1929198 , 755.68207248],
[755.10969541, 756.38178297, 761.18044623, ..., 762.22995695,
756.69732237, 756.35418424],
[753.32406793, 752.80673891, 752.03157068, ..., 756.22896884,
755.6715562 , 756.78221912]])array([[0.00089158, 0.00087404, 0.0008765 , ..., 0.00088058, 0.00044662,
0.00087892],
[0.00087075, 0.00088254, 0.00088179, ..., 0.00086737, 0.00170292,
0.0008685 ],
[0.00089083, 0.0008845 , 0.00088362, ..., 0.00087608, 0.00044392,
0.00086608],
[0.00096813, 0.00086833, 0.00168521, ..., 0.00089063, 0.00129646,
0.00166637]])array([[3, 3, 3, ..., 3, 2, 3],
[3, 3, 3, ..., 3, 4, 3],
[3, 3, 3, ..., 3, 2, 3],
[3, 3, 4, ..., 3, 4, 4]])array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]])array([[0.41500359, 0.41500359, 0.41500359, ..., 0.41500359, 0.41500359,
0.41500359],
[0.42993239, 0.42993239, 0.42993239, ..., 0.42993239, 0.42993239,
0.42993239],
[0.55799239, 0.55799239, 0.55799239, ..., 0.55799239, 0.55799239,
0.55799239],
[0.39177214, 0.39177214, 0.39177214, ..., 0.39177214, 0.39177214,
0.39177214]])array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]])array([[53460.36288858, 53460.36385433, 53460.36478871, ...,
53461.40632654, 53461.40726571, 53461.40777038],
[53460.37427417, 53460.37521658, 53460.37615904, ...,
53461.47633108, 53461.47725754, 53461.479019 ],
[53460.42243175, 53460.42339404, 53460.42433825, ...,
53461.50554 , 53461.50647421, 53461.50697612],
[53460.39759808, 53460.39864992, 53460.39957983, ...,
53461.51435183, 53461.51530904, 53461.51666683]])array([[0.75118358, 0.54917103, 0.90600545, ..., 0.64640509, 0.24651666,
0.84703996],
[1. , 1. , 0.18921988, ..., 0.81098596, 0.34035628,
0.99696537],
[0.56560841, 0.95846614, 0.85780695, ..., 0.50105863, 0.99947793,
1. ],
[0.9683953 , 0.61509461, 0.99814479, ..., 0.9439049 , 0.99189283,
0.82635188]])array([[ 6, 4, -6, ..., 4, 2, 5],
[ 4, 7, 0, ..., -6, 9, -2],
[ 3, 5, 5, ..., -4, -2, 4],
[-3, -4, -3, ..., 4, -1, 9]])PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))
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))<xarray.Dataset>
Dimensions: (count_obs: 250)
Coordinates:
* count_obs (count_obs) int64 0 1 2 3 4 5 6 7 ... 243 244 245 246 247 248 249
Data variables:
count (count_obs) int64 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: 2023-12-05T18:56:31.600135
arviz_version: 0.16.1
inference_library: pymc
inference_library_version: 5.8.1
modeling_interface: bambi
modeling_interface_version: 0.13.0.dev0array([ 0, 1, 2, ..., 247, 248, 249])
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])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='count_obs', length=250))draw_1 = select_draws(idata_grid, grid, {"livebait": 0}, "count_mean")
draw_2 = select_draws(idata_grid, grid, {"livebait": 1}, "count_mean")
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")["count_mean"].values,
"hdi_high": comparison_hdi.sel(hdi="higher")["count_mean"].values,
}
)
comparison_df.head(10)| mean | hdi_low | hdi_high | |
|---|---|---|---|
| 0 | 0.214363 | 0.144309 | 0.287735 |
| 1 | 0.053678 | 0.029533 | 0.077615 |
| 2 | 0.013558 | 0.006332 | 0.021971 |
| 3 | 0.003454 | 0.001132 | 0.006040 |
| 4 | 0.512709 | 0.369741 | 0.661034 |
| 5 | 0.128316 | 0.077068 | 0.181741 |
| 6 | 0.032392 | 0.015553 | 0.050690 |
| 7 | 0.008247 | 0.003047 | 0.014382 |
| 8 | 1.228708 | 0.913514 | 1.533121 |
| 9 | 0.307342 | 0.192380 | 0.426808 |
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.214363 | 0.144309 | 0.287735 |
| 1 | livebait | diff | (0, 1) | 0 | 1 | 2 | 0.053678 | 0.029533 | 0.077615 |
| 2 | livebait | diff | (0, 1) | 0 | 1 | 3 | 0.013558 | 0.006332 | 0.021971 |
| 3 | livebait | diff | (0, 1) | 0 | 1 | 4 | 0.003454 | 0.001132 | 0.006040 |
| 4 | livebait | diff | (0, 1) | 0 | 2 | 1 | 0.512709 | 0.369741 | 0.661034 |
| 5 | livebait | diff | (0, 1) | 0 | 2 | 2 | 0.128316 | 0.077068 | 0.181741 |
| 6 | livebait | diff | (0, 1) | 0 | 2 | 3 | 0.032392 | 0.015553 | 0.050690 |
| 7 | livebait | diff | (0, 1) | 0 | 2 | 4 | 0.008247 | 0.003047 | 0.014382 |
| 8 | livebait | diff | (0, 1) | 0 | 3 | 1 | 1.228708 | 0.913514 | 1.533121 |
| 9 | livebait | diff | (0, 1) | 0 | 3 | 2 | 0.307342 | 0.192380 | 0.426808 |
Albeit the other information in the summary_df, the columns estimate, lower_3.0%, upper_97.0% are identical.
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.
| term | estimate_type | value | camper | persons | child | estimate | lower_3.0% | upper_97.0% | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | livebait | diff | (0, 1) | 0 | 1 | 1 | 0.214363 | 0.144309 | 0.287735 |
| 1 | livebait | diff | (0, 1) | 0 | 1 | 2 | 0.053678 | 0.029533 | 0.077615 |
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, "count_mean")
draws_11 = select_draws(idata_grid, grid, cond_11, "count_mean")
(draws_11 - draws_10).mean(("chain", "draw")).item()0.2143627093182434Next, we need to compute the comparison for row 1.
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.16068445232635978To verify this is correct, we can check by performing the cross-comparison directly on the summary_df.
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.