Implementing instrumental variable and backdoor criteria plus propensity and regression methods
Notes
- Causal Bayesian Networkx: here
dowhy > causal_model.py CausalModel > causal_graph.py CausalGraph
%matplotlib inline
%load_ext autoreload
%autoreload 2
plt.style.use('bmh')
treatments = ['V0', 'V1']
outcome = 'Y'
common_causes = ['W0']
effect_modifiers = ['X0']
instruments = []
observed_nodes = treatments + [outcome] + instruments + effect_modifiers + common_causes
add_unobserved_confounder = True
missing_nodes_as_confounders = True
cg_ref = dw.causal_graph.CausalGraph(treatments,
[outcome], graph=None,
common_cause_names=common_causes,
instrument_names=instruments,
effect_modifier_names=effect_modifiers,
observed_node_names=observed_nodes)
cg_ref._graph.nodes(data=True)
cg = CausalGraph(treatments=treatments, outcome=outcome, common_causes=common_causes,
effect_modifiers=effect_modifiers, observed_nodes=observed_nodes,
missing_nodes_as_confounders=missing_nodes_as_confounders,
add_unobserved_confounder=add_unobserved_confounder)
cg.g.nodes['U']['observed']
cg.view_graph()
cg.get_ancestors('V0')
g_cut = cg.cut_edges([('U','Y'), ('W0', 'V1')])
show_graph(g_cut)
cg.get_causes(['V0'])
cg.get_instruments(treatments, outcome)
cg.get_effect_modifiers(treatments, [outcome])
treatments = ['V0',] # 'V1']
outcome = 'Y'
common_causes = ['W0']
effect_modifiers = ['X0']
instruments = []
observed_nodes = treatments + [outcome] + instruments + effect_modifiers
add_unobserved_confounder = True
missing_nodes_as_confounders = True
cg_kwargs = dict(
missing_nodes_as_confounders=missing_nodes_as_confounders,
add_unobserved_confounder=add_unobserved_confounder,
observed_nodes=observed_nodes
)
cm = CausalModel(treatments=treatments, outcome=outcome, common_causes=common_causes,
effect_modifiers=effect_modifiers,
causal_graph_kwargs=cg_kwargs)
estimands = cm.identify_effect(); estimands
Regression estimators based on sklearn regression classes
isinstance(linear_model.LinearRegression(), sklearn.base.RegressorMixin)
Sanity checking on a quadratic polynomial toy dataset
X = np.linspace(-1, 1, 200)
X = np.array([X**2, X, np.ones(len(X))*.5]).T
w = np.array([2, 0, .5])
y = X @ w
fig, ax = plt.subplots(figsize=(8,4), constrained_layout=True)
ax.scatter(X[:,-2], y)
ax.set(xlabel='x', ylabel='y', title='dataset')
plt.show()
regression_model = linear_model.LinearRegression()
estimator = RegressionEstimator(regression_model)
estimator.fit(X, y, ix=0, ix_confounders=[1])
ate = estimator.estimate_effect(X=X, treatment=1, control=0)
print(f'ate = {ate:.3f} coefs {estimator.m.coef_}')
Classification estimator
propensity score: common causes -> prediction of treatment (class) -> grouping by score to select pairs of most similar treatment and control group samples to compute the difference in outcome
grouping is done using some nearest neighbour search:
- ATC if nearest neighbor is set up with the treated group and for each control group sample a match is looked up and then the difference of the outcome is computed
- ATT if nearest neighbor is set up with the control group and for each treated group sample a match is looked up and then the difference of the outcome is computed
TODO: test PropensityScoreMatcher on data generated using bcg.basics classes
n = 200
x_treatment = np.random.choice([True, False], p=[.5, .5], size=n)
x_common_causes = np.array([
[np.random.normal(loc=v, scale=.1) for v in x_treatment],
[np.random.normal(loc=10-v, scale=.1) for v in x_treatment],
])
y_outcome = np.array([np.random.normal(loc=v, scale=.1) for v in x_treatment])
fig, axs = plt.subplots(figsize=(8,6), nrows=4, constrained_layout=True)
axs[0].hist(x_treatment.astype(float))
axs[0].set(xlabel='treatment')
axs[1].hist(x_common_causes[0])
axs[1].set(xlabel='cc0')
axs[2].hist(x_common_causes[1])
axs[2].set(xlabel='cc1')
axs[3].hist(y_outcome)
axs[3].set(xlabel='outcome')
plt.show()
fig, ax = plt.subplots(figsize=(8,4), constrained_layout=True)
ax.scatter(x_treatment, y_outcome)
ax.set(xlabel='treatment', ylabel='outcome', title='dataset')
plt.show()
X, y = np.concatenate((x_treatment[:,None], x_common_causes.T), axis=1), y_outcome
X.shape, y.shape
class PropensityScoreMatcher:
def __init__(self, propensity_model:sklearn.base.ClassifierMixin):
assert isinstance(propensity_model, sklearn.base.ClassifierMixin)
self.pm = propensity_model
def fit(self, X:np.ndarray, y:np.ndarray, ix:int, ix_confounders:List[int], reset:bool=True):
'''building the classifier model & nearest neigbhor search thingy
ix: needs to point to a binary variable
'''
if not isinstance(ix_confounders, list):
ix_confounders = list(ix_confounders)
self.ix = ix
self.ix_confounders = ix_confounders
_ix = [ix] + ix_confounders
self._ix = _ix
if reset:
self.pm.fit(X[:, self.ix_confounders], X[:,self.ix])
def estimate_effect(self, X:np.ndarray, treatment:Union[int, bool], control:Union[int, bool],
y:np.ndarray=None, kind:str='ate'):
assert y is not None, 'Cannot be None. That\'s just the default to have consistent method parameters.'
assert kind in ['ate', 'att', 'atc']
propensity_score = self.pm.predict(X[:, self.ix_confounders])
ix_treat, ix_control = X[:,self.ix] == treatment, X[:,self.ix] == control
X_treat, X_cont = X[ix_treat,:], X[ix_control,:]
y_treat, y_cont = y[ix_treat], y[ix_control]
searcher = neighbors.NearestNeighbors(n_neighbors=1)
def get_att():
searcher.fit(propensity_score[ix_control][:,None])
distances, indices = searcher.kneighbors(propensity_score[ix_treat][:,None])
att = 0
n_treat = ix_treat.sum()
for i in range(n_treat):
out_treat = y_treat[i]
out_cont = y_cont[indices[i][0]]
att += out_treat - out_cont
return att / n_treat
def get_atc():
searcher.fit(propensity_score[ix_treat][:,None])
distances, indices = searcher.kneighbors(propensity_score[ix_control][:,None])
atc = 0
n_cont = ix_control.sum()
for i in range(n_cont):
out_treat = y_treat[indices[i][0]]
out_cont = y_cont[i]
atc += out_treat - out_cont
return atc / n_cont
def get_ate():
n_treat = ix_treat.sum()
n_cont = ix_control.sum()
att = get_att()
atc = get_atc()
return (att*n_treat + atc*n_cont) / (n_treat + n_cont)
if kind == 'ate':
return get_ate()
elif kind == 'att':
return get_att()
elif kind == 'atc':
return get_atc()
else:
raise NotImplementedError
propensity_model = linear_model.LogisticRegression(solver='lbfgs')
estimator = PropensityScoreMatcher(propensity_model)
estimator.fit(X, y, ix=0, ix_confounders=[1, 2])
ate = estimator.estimate_effect(X=X, treatment=True, control=False, y=y)
print(f'ate = {ate:.3f}')
Generating data for the graphical model using bcg.basics functions
outcome_is_binary = True
treatment_is_binary = True
n = 333
n_common_causes = len(common_causes)
n_instruments = len(instruments)
n_eff_mods = len(effect_modifiers)
n_treatments = len(treatments)
beta = 1 # max random value
cc = CommonCauses.get_obs(n, n_common_causes)
ins = Instruments.get_obs(n, n_instruments)
em = EffectModifiers.get_obs(n, n_eff_mods)
treat = Treatments.get_obs(n, n_treatments, cc, ins, beta, treatment_is_binary=treatment_is_binary)
out = Outcomes.get_obs(treat, cc, em, outcome_is_binary=outcome_is_binary)
obs = pd.concat((treat.obs, cc.obs, em.obs, ins.obs, out.obs), axis=1)
X, y, not_target = get_Xy(obs, target=outcome)
obs.head(), obs.tail()
not_target.index('V0')
Adding effect estimate functionality to CausalModel
Changing the implementation of get_Xy, incorporating products with effect modifiers, based on lns 59-71 in causal_estimators/linear_regression_estimator.py with the new argumentfeature_product_groups. The variable is supposed to consist of two lists, each containing features in obs, of which all products will be computed.
get_Xy_with_products(obs, target=outcome, feature_product_groups=[treatments, effect_modifiers])
causal_method = 'backdoor'
control_value = 0
treatment_name = 'v0'
treatment_value = 2
effect_modifiers = effect_modifiers
target_unit = 'ate'
# model = linear_model.LinearRegression()
# model = linear_model.LogisticRegression()
model = None
supervised_type_is_regression = False
cm.estimate_effect(estimands, control_value, treatment_name, treatment_value,
obs, outcome=outcome, causal_method=causal_method, model=model,
target_unit=target_unit, effect_modifiers=effect_modifiers,
supervised_type_is_regression=supervised_type_is_regression)
a = np.linspace(1, 4, 5)
b = a[:, np.newaxis]; b