HierarchicalLinearRegression#

class causalpy.pymc_models.HierarchicalLinearRegression[source]#

Custom PyMC model for Hierarchical linear regression with fixed and group-level random effects.

Defines the PyMC model (centered parameterization)

\[\begin{split}\beta_{\text{fixed}} &\sim \mathrm{Normal}(0, 10) \\ \sigma_{\text{random}} &\sim \mathrm{HalfNormal}(1) \\ \beta_{\text{random}, g} &\sim \mathrm{Normal}(0, \sigma_{\text{random}}) \\ \mu &= X \cdot \beta_{\text{fixed}}^{\top} + \sum\left(\beta_{\text{random}}[\mathrm{group\_idx}] \odot Z, \mathrm{axis}=2\right) \\ \sigma_{\text{fixed}} &\sim \mathrm{HalfNormal}(1) \\ y &\sim \mathrm{Normal}(\mu, \sigma_{\text{fixed}})\end{split}\]

Set priors with keys beta_fixed, sigma_random, and sigma_fixed.

For the non-centered parameterization:

\[\begin{split}\beta_{\text{group}} &\sim \mathrm{Normal}(0, 1) \\ \beta_{\text{random}} &= \beta_{\text{group}} \odot \sigma_{\text{random}}\end{split}\]

Set priors with keys beta_fixed, sigma_random, sigma_fixed, and beta_group.

Example

>>> import numpy as np
>>> import xarray as xr
>>> import pymc as pm
>>> rng = np.random.default_rng(42)
>>> n_groups, n_obs = 4, 40
>>> group_idx = np.repeat(np.arange(n_groups), n_obs // n_groups).astype(np.int32)
>>> x1 = rng.normal(size=n_obs)
>>> X = xr.DataArray(
...     np.column_stack([np.ones(n_obs), x1]),
...     dims=["obs_ind", "coeffs"],
...     coords={"obs_ind": np.arange(n_obs), "coeffs": ["1", "x1"]},
... )
>>> Z = xr.DataArray(
...     np.column_stack([np.ones(n_obs), x1]),
...     dims=["obs_ind", "random_coeffs"],
...     coords={
...         "obs_ind": np.arange(n_obs),
...         "random_coeffs": ["1|group", "x1|group"],
...     },
... )
>>> y = xr.DataArray(
...     rng.normal(size=(n_obs, 1)),
...     dims=["obs_ind", "treated_units"],
...     coords={"obs_ind": np.arange(n_obs), "treated_units": ["unit_0"]},
... )
>>> coords = {
...     "obs_ind": np.arange(n_obs),
...     "coeffs": ["1", "x1"],
...     "random_coeffs": ["1|group", "x1|group"],
...     "treated_units": ["unit_0"],
...     "groups": np.arange(n_groups),
... }
>>> priors = {
...     "beta_fixed": Prior("Normal", mu=0, sigma=10, dims=["treated_units", "coeffs"]),
...     "sigma_fixed": Prior("HalfNormal", sigma=1, dims=["treated_units"]),
...     "sigma_random": Prior("HalfNormal", sigma=1, dims=["treated_units", "random_coeffs"]),
...     "beta_group": Prior("Normal", mu=0, sigma=1, dims=["groups", "treated_units", "random_coeffs"]),
... }
>>> model = HierarchicalLinearRegression(
...     sample_kwargs={"draws": 10, "tune": 10, "chains": 1, "progressbar": False},
...     priors=priors,
... )
>>> model.build_model(
...     X=X,
...     Z=Z,
...     y=y,
...     group_idx=group_idx,
...     coords=coords,
...     non_centered=True,
... )
>>> with model:
...     idata = pm.sample(**model.sample_kwargs)
>>> "mu" in idata.posterior
True

Methods

`HierarchicalLinearRegression.add_coord`(name)	Register a dimension coordinate with the model.
`HierarchicalLinearRegression.add_coords`(...)	Vectorized version of `Model.add_coord`.
`HierarchicalLinearRegression.add_named_variable`(var)	Add a random graph variable to the named variables of the model.
`HierarchicalLinearRegression.build_model`(X, ...)	Defines the PyMC model.
`HierarchicalLinearRegression.calculate_cumulative_impact`(impact)
`HierarchicalLinearRegression.calculate_impact`(...)	Calculate the causal impact as the difference between observed and predicted values.
`HierarchicalLinearRegression.check_start_vals`(...)	Check that the logp is defined and finite at the starting point.
`HierarchicalLinearRegression.compile_d2logp`([...])	Compiled log probability density hessian function.
`HierarchicalLinearRegression.compile_dlogp`([...])	Compiled log probability density gradient function.
`HierarchicalLinearRegression.compile_fn`(outs, *)
`HierarchicalLinearRegression.compile_logp`([...])	Compiled log probability density function.
`HierarchicalLinearRegression.copy`()	Clone the model.
`HierarchicalLinearRegression.create_value_var`(...)	Create a `TensorVariable` that will be used as the random variable's "value" in log-likelihood graphs.
`HierarchicalLinearRegression.d2logp`([vars, ...])	Hessian of the models log-probability w.r.t.
`HierarchicalLinearRegression.debug`([point, ...])	Debug model function at point.
`HierarchicalLinearRegression.dlogp`([vars, ...])	Gradient of the models log-probability w.r.t.
`HierarchicalLinearRegression.eval_rv_shapes`()	Evaluate shapes of untransformed AND transformed free variables.
`HierarchicalLinearRegression.fit`(X, Z, y, ...)	Draw posterior and predictive samples for hierarchical linear regression.
`HierarchicalLinearRegression.get_context`([...])
`HierarchicalLinearRegression.initial_point`([...])	Compute the initial point of the model.
`HierarchicalLinearRegression.logp`([vars, ...])	Elemwise log-probability of the model.
`HierarchicalLinearRegression.logp_dlogp_function`([...])	Compile a PyTensor function that computes logp and gradient.
`HierarchicalLinearRegression.make_obs_var`(...)	Create a TensorVariable for an observed random variable.
`HierarchicalLinearRegression.name_for`(name)	Check if name has prefix and adds if needed.
`HierarchicalLinearRegression.name_of`(name)	Check if name has prefix and deletes if needed.
`HierarchicalLinearRegression.point_logps`([...])	Compute the log probability of point for all random variables in the model.
`HierarchicalLinearRegression.predict`(X[, ...])	Predict data given input data X
`HierarchicalLinearRegression.print_coefficients`(labels)	Print the model coefficients with their labels.
`HierarchicalLinearRegression.priors_from_data`(X, y)	Generate priors dynamically based on the input data.
`HierarchicalLinearRegression.profile`(outs, *)	Compile and profile a PyTensor function which returns `outs` and takes values of model vars as a dict as an argument.
`HierarchicalLinearRegression.register_data_var`(data)	Register a data variable with the model.
`HierarchicalLinearRegression.register_rv`(...)	Register an (un)observed random variable with the model.
`HierarchicalLinearRegression.replace_rvs_by_values`(...)	Clone and replace random variables in graphs with their value variables.
`HierarchicalLinearRegression.score`(X, y[, ...])	Score the Bayesian \(R^2\) given inputs `X` and outputs `y`.
`HierarchicalLinearRegression.set_data`(name, ...)	Change the values of a data variable in the model.
`HierarchicalLinearRegression.set_dim`(name, ...)	Update a mutable dimension.
`HierarchicalLinearRegression.set_initval`(...)	Set an initial value (strategy) for a random variable.
`HierarchicalLinearRegression.shape_from_dims`(dims)
`HierarchicalLinearRegression.to_graphviz`(*)	Produce a graphviz Digraph from a PyMC model.

Attributes

`basic_RVs`	List of random variables the model is defined in terms of.
`continuous_value_vars`	All the continuous value variables in the model.
`coords`	Coordinate values for model dimensions.
`datalogp`	PyTensor scalar of log-probability of the observed variables and potential terms.
`default_priors`
`dim_lengths`	The symbolic lengths of dimensions in the model.
`discrete_value_vars`	All the discrete value variables in the model.
`isroot`
`observedlogp`	PyTensor scalar of log-probability of the observed variables.
`parent`
`potentiallogp`	PyTensor scalar of log-probability of the Potential terms.
`prefix`
`root`
`unobserved_RVs`	List of all random variables, including deterministic ones.
`unobserved_value_vars`	List of all random variables (including untransformed projections), as well as deterministics used as inputs and outputs of the model's log-likelihood graph.
`value_vars`	List of unobserved random variables used as inputs to the model's log-likelihood (which excludes deterministics).
`varlogp`	PyTensor scalar of log-probability of the unobserved random variables (excluding deterministic).
`varlogp_nojac`	PyTensor scalar of log-probability of the unobserved random variables (excluding deterministic) without jacobian term.

__init__(sample_kwargs=None, priors=None)#

Parameters:

sample_kwargs (dict[str, Any] | None) – Dictionary of kwargs that get unpacked and passed to the pymc.sample() function. Defaults to an empty dictionary if None.
priors (dict[str, Any] | None) – Dictionary of priors for the model. Defaults to None, in which case default priors are used.

Return type:

None

classmethod __new__(*args, **kwargs)#