InterventionTimeEstimator#

class causalpy.pymc_models.InterventionTimeEstimator[source]#

Custom PyMC model to estimate the time an intervention took place.

This model implements three types of changepoints: level shift, trend change, and impulse response. While the underlying mathematical framework could theoretically be applied to other changepoint detection problems, it has been specifically designed and tested for use in interrupted time series causal inference settings.

In words: - beta represents the regression coefficients for the baseline signal μ. - tau is the changepoint time, and w is a sigmoid function to smooth the transition. - level and trend model a linear shift after the changepoint. - A and lambda define an impulse response at the changepoint. - mu_in combines the level, trend, and impulse contributions. - mu_ts is the total mean, including baseline and intervention. - sigma is the observation noise. - Finally, y is drawn from a Normal distribution with mean mu_ts and standard deviation sigma.

\[\begin{split}\beta &\sim \mathrm{Normal}(0, 5) \\ \mu &= \beta \cdot X\\ \\ \tau &\sim \mathrm{Uniform}(\text{lower_bound}, \text{upper_bound}) \\ w &= sigmoid(t-\tau) \\ \\ \text{level} &\sim \mathrm{Normal}(0, 5) \\ \text{trend} &\sim \mathrm{Normal}(0, 0.5) \\ A &\sim \mathrm{Normal}(0, 5) \\ \lambda &\sim \mathrm{HalfNormal}(0, 5) \\ \text{impulse} &= A \cdot exp(-\lambda \cdot |t-\tau|) \\ \mu_{in} &= \text{level} + \text{trend} \cdot (t-\tau) + \text{impulse}\\ \\ \sigma &\sim \mathrm{HalfNormal}(0, 1) \\ \mu_{ts} &= \mu + \mu_{in} \\ \\ y &\sim \mathrm{Normal}(\mu_{ts}, \sigma)\end{split}\]

Example

>>> import causalpy as cp
>>> import numpy as np
>>> from patsy import build_design_matrices, dmatrices
>>> from causalpy.pymc_models import InterventionTimeEstimator as ITE
>>> data = cp.load_data("its")
>>> formula="y ~ 1 + t + C(month)"
>>> y, X = dmatrices(formula, data)
>>> outcome_variable_name = y.design_info.column_names[0]
>>> labels = X.design_info.column_names
>>> _y, _X = np.asarray(y), np.asarray(X)
>>> _X = xr.DataArray(
...     _X,
...     dims=["obs_ind", "coeffs"],
...     coords={
...         "obs_ind": data.index,
...         "coeffs": labels,
...         },
... )
>>> _y = xr.DataArray(
...     _y,
...     dims=["obs_ind", "treated_units"],
...     coords={
...         "obs_ind": data.index,
...         "treated_units": ["unit_0"]
...         },
... )
>>> COORDS = {
...     "coeffs": labels,
...     "obs_ind": np.arange(X.shape[0]),
...     "treated_units": ["unit_0"],
... }
>>> model = ITE(treatment_effect_type="level", sample_kwargs={"draws" : 10, "tune":10, "progressbar":False})
>>> model.set_time_range(None, data)
>>> model.fit(X=_X, y=_y, coords=COORDS)
Inference ...

Methods

`InterventionTimeEstimator.__init__`(...[, ...])	Initializes the InterventionTimeEstimator model.
`InterventionTimeEstimator.add_coord`(name[, ...])	Register a dimension coordinate with the model.
`InterventionTimeEstimator.add_coords`(coords, *)	Vectorized version of `Model.add_coord`.
`InterventionTimeEstimator.add_named_variable`(var)	Add a random graph variable to the named variables of the model.
`InterventionTimeEstimator.build_model`(X, y, ...)	Defines the PyMC model
`InterventionTimeEstimator.calculate_cumulative_impact`(impact)
`InterventionTimeEstimator.calculate_impact`(...)
`InterventionTimeEstimator.check_start_vals`(...)	Check that the logp is defined and finite at the starting point.
`InterventionTimeEstimator.compile_d2logp`([...])	Compiled log probability density hessian function.
`InterventionTimeEstimator.compile_dlogp`([...])	Compiled log probability density gradient function.
`InterventionTimeEstimator.compile_fn`(outs, *)	Compiles a PyTensor function.
`InterventionTimeEstimator.compile_logp`([...])	Compiled log probability density function.
`InterventionTimeEstimator.copy`()	Clone the model.
`InterventionTimeEstimator.create_value_var`(...)	Create a `TensorVariable` that will be used as the random variable's "value" in log-likelihood graphs.
`InterventionTimeEstimator.d2logp`([vars, ...])	Hessian of the models log-probability w.r.t.
`InterventionTimeEstimator.debug`([point, fn, ...])	Debug model function at point.
`InterventionTimeEstimator.dlogp`([vars, jacobian])	Gradient of the models log-probability w.r.t.
`InterventionTimeEstimator.eval_rv_shapes`()	Evaluate shapes of untransformed AND transformed free variables.
`InterventionTimeEstimator.fit`(X, y[, coords])	Draw samples from posterior, prior predictive, and posterior predictive distributions, placing them in the model's idata attribute.
`InterventionTimeEstimator.get_context`([...])
`InterventionTimeEstimator.initial_point`([...])	Compute the initial point of the model.
`InterventionTimeEstimator.logp`([vars, ...])	Elemwise log-probability of the model.
`InterventionTimeEstimator.logp_dlogp_function`([...])	Compile a PyTensor function that computes logp and gradient.
`InterventionTimeEstimator.make_obs_var`(...)	Create a TensorVariable for an observed random variable.
`InterventionTimeEstimator.name_for`(name)	Check if name has prefix and adds if needed.
`InterventionTimeEstimator.name_of`(name)	Check if name has prefix and deletes if needed.
`InterventionTimeEstimator.point_logps`([...])	Compute the log probability of point for all random variables in the model.
`InterventionTimeEstimator.predict`(X)	Predict data given input data X
`InterventionTimeEstimator.print_coefficients`(labels)
`InterventionTimeEstimator.profile`(outs, *[, ...])	Compile and profile a PyTensor function which returns `outs` and takes values of model vars as a dict as an argument.
`InterventionTimeEstimator.register_data_var`(data)	Register a data variable with the model.
`InterventionTimeEstimator.register_rv`(...[, ...])	Register an (un)observed random variable with the model.
`InterventionTimeEstimator.replace_rvs_by_values`(...)	Clone and replace random variables in graphs with their value variables.
`InterventionTimeEstimator.score`(X, y)	Score the Bayesian \(R^2\) given inputs `X` and outputs `y`.
`InterventionTimeEstimator.set_data`(name, values)	Change the values of a data variable in the model.
`InterventionTimeEstimator.set_dim`(name, ...)	Update a mutable dimension.
`InterventionTimeEstimator.set_initval`(...)	Set an initial value (strategy) for a random variable.
`InterventionTimeEstimator.set_time_range`(...)	Set time_range.
`InterventionTimeEstimator.shape_from_dims`(dims)
`InterventionTimeEstimator.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.
`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__(treatment_effect_type, treatment_effect_param=None, sample_kwargs=None)[source]#

Initializes the InterventionTimeEstimator model.

Parameters:

treatment_effect_type (str | list[str]) –
Optional dictionary that specifies prior parameters for the intervention effects. Expected keys are:
- ”level”: [mu, sigma]
- ”trend”: [mu, sigma]
- ”impulse”: [mu, sigma1, sigma2]
If a key is missing, the corresponding effect is ignored. If the associated list is incomplete, default values will be used.
sample_kwargs – Optional dictionary of arguments passed to pm.sample().

classmethod __new__(*args, **kwargs)#