Causal Effect Estimation in datasets

This notebook serves as a demonstration of the CausalEffectEstimation module for estimating causal effects on both generated and real datasets.

import numpy as np
import pandas as pd

import pyagrum as gum
import pyagrum.lib.discretizer as disc
import pyagrum.lib.notebook as gnb

import pyagrum.explain as gexpl

import pyagrum.causalEffectEstimation as gcee

1. ACE estimations from generated RCT data

To estimate the Average Causal Effects (ACE) from data generated under conditions resembling a Randomized Controlled Trial (RCT), we define two generative models.

1.1 Dataset

Consider two generative models:

A linear generative model described by the equation: $Y = 3X_1 + 2X_2 -2X_3 -0.8X_4 + T(2X_1 + 5X_3 +3X_4)$
And a non-linear generative model described by the equation: $Y = 3X_1 + 2X_2^2 -2X_3 -0.8X_4 +10T$

Where $(X_1,X_2,X_3,X_4) \sim \mathcal{N}_4((1,1,1,1), I_4)$ , $T \sim \mathcal{Ber}(1/2)$ and $(X_1,X_2,X_3,X_4,T)$ are jointly independent in both of the models.

Data from the models can be obtatined by the functions given below.

def linear_simulation(n: int = 100000, sigma: float = 1) -> pd.DataFrame:
  """
  Returns n observations from the linear model with normally distributed
  noise with expected value 0 and standard deviation sigma.
  """

  X1 = np.random.normal(1, 1, n)
  X2 = np.random.normal(1, 1, n)
  X3 = np.random.normal(1, 1, n)
  X4 = np.random.normal(1, 1, n)
  epsilon = np.random.normal(0, sigma, n)
  T = np.random.binomial(1, 0.5, n)
  Y = 3 * X1 + 2 * X2 - 2 * X3 - 0.8 * X4 + T * (2 * X1 + 5 * X3 + 3 * X4) + epsilon
  d = np.array([T, X1, X2, X3, X4, Y])
  df_data = pd.DataFrame(data=d.T, columns=["T", "X1", "X2", "X3", "X4", "Y"])
  df_data["T"] = df_data["T"].astype(int)

  return df_data


def non_linear_simulation(n: int = 100000, sigma: float = 1) -> pd.DataFrame:
  """
  Returns n observations from the non-linear model with normally distributed
  noise with expected value 0 and standard deviation sigma.
  """

  X1 = np.random.normal(1, 1, n)
  X2 = np.random.normal(1, 1, n)
  X3 = np.random.normal(1, 1, n)
  X4 = np.random.normal(1, 1, n)
  epsilon = np.random.normal(0, sigma, n)
  T = np.random.binomial(1, 0.5, n)
  Y = 3 * X1 + 2 * X2**2 - 2 * X3 - 0.8 * X4 + 10 * T + epsilon
  d = np.array([T, X1, X2, X3, X4, Y])
  df_data = pd.DataFrame(data=d.T, columns=["T", "X1", "X2", "X3", "X4", "Y"])
  df_data["T"] = df_data["T"].astype(int)

  return df_data

Furthermore, the expected values of $Y(0)$ and $Y(1)$ can be explicitly calculated, providing us the theoretical ACE of $\tau = 10$ which will serve as a point of reference for the estimations.

We will explore how the CausalEffectEstimation module can estimate the causal effect of $T$ on $Y$ in both of the generated datasets.

1.2 Setup

First, genereate the data using the previously defined functions, and specify the causal graph of the variables. A single graph will be applicable to both datasets, as they share the same variables and causal structure.

linear_data = linear_simulation()
non_linear_data = non_linear_simulation()

bn = gum.fastBN("Y<-T; Y<-X1; Y<-X2; Y<-X3; Y<-X4")
causal_model = gum.CausalModel(bn)

gnb.showCausalModel(causal_model, size="10")

svg

If the CausalBNEstimator estimation method is not employed for estimation process, the model generated using gum.fastBN will suffice, as only the graph structure of the Causal Bayesian Network will be used for causal identification.

However, if the Causal Bayesian Network estimator is utilized, it will be necessary to provide a gum.CausalModel object with appropriate discretization, as the Conditional Probability Tables of the model will be used for estimation. Here we use the discretizer module to perform this task, the arcs are added manually.

Selecting an optimal discretization is crucial: a coarse discretization may lead to poor estimation due to its inability to capture fine variations in the distribution, while an overly fine discretization may result in too many parameters, making it difficult for the parameter learning algorithm to accurately estimate the distribution. Therefore, the discretization should strike a balance, being neither too coarse nor too fine.

discretizer = disc.Discretizer(defaultDiscretizationMethod="uniform", defaultNumberOfBins=20)
disc_bn = discretizer.discretizedTemplate(linear_data)

disc_bn.beginTopologyTransformation()
for node in {"T", "X1", "X2", "X3", "X4"}:
  disc_bn.addArc(node, "Y")
disc_bn.endTopologyTransformation()

disc_causal_model = gum.CausalModel(disc_bn)

for id, _ in disc_bn:
  print(disc_bn.variable(id))

T:Range([0,1])
X1:Discretized(<(-3.5643869716843213;-3.109312457102246[,[-3.109312457102246;-2.6542379425201705[,[-2.6542379425201705;-2.199163427938095[,[-2.199163427938095;-1.7440889133560198[,[-1.7440889133560198;-1.2890143987739444[,[-1.2890143987739444;-0.833939884191869[,[-0.833939884191869;-0.37886536960979367[,[-0.37886536960979367;0.0762091449722817[,[0.0762091449722817;0.5312836595543571[,[0.5312836595543571;0.9863581741364325[,[0.9863581741364325;1.4414326887185078[,[1.4414326887185078;1.8965072033005832[,[1.8965072033005832;2.3515817178826586[,[2.3515817178826586;2.806656232464734[,[2.806656232464734;3.2617307470468093[,[3.2617307470468093;3.7168052616288847[,[3.7168052616288847;4.17187977621096[,[4.17187977621096;4.626954290793035[,[4.626954290793035;5.08202880537511[,[5.08202880537511;5.537103319957185)>)
X2:Discretized(<(-4.441687990260335;-3.9481321348423655[,[-3.9481321348423655;-3.4545762794243955[,[-3.4545762794243955;-2.961020424006426[,[-2.961020424006426;-2.467464568588456[,[-2.467464568588456;-1.9739087131704864[,[-1.9739087131704864;-1.480352857752517[,[-1.480352857752517;-0.9867970023345469[,[-0.9867970023345469;-0.4932411469165774[,[-0.4932411469165774;0.0003147085013921469[,[0.0003147085013921469;0.4938705639193621[,[0.4938705639193621;0.9874264193373321[,[0.9874264193373321;1.4809822747553012[,[1.4809822747553012;1.9745381301732712[,[1.9745381301732712;2.468093985591241[,[2.468093985591241;2.9616498410092102[,[2.9616498410092102;3.4552056964271802[,[3.4552056964271802;3.9487615518451493[,[3.9487615518451493;4.442317407263119[,[4.442317407263119;4.935873262681089[,[4.935873262681089;5.429429118099059)>)
X3:Discretized(<(-3.29718341748295;-2.8557020243544247[,[-2.8557020243544247;-2.4142206312259[,[-2.4142206312259;-1.9727392380973745[,[-1.9727392380973745;-1.5312578449688494[,[-1.5312578449688494;-1.0897764518403243[,[-1.0897764518403243;-0.648295058711799[,[-0.648295058711799;-0.2068136655832742[,[-0.2068136655832742;0.23466772754525111[,[0.23466772754525111;0.6761491206737764[,[0.6761491206737764;1.1176305138023013[,[1.1176305138023013;1.5591119069308261[,[1.5591119069308261;2.000593300059352[,[2.000593300059352;2.4420746931878767[,[2.4420746931878767;2.8835560863164016[,[2.8835560863164016;3.3250374794449273[,[3.3250374794449273;3.766518872573452[,[3.766518872573452;4.208000265701977[,[4.208000265701977;4.649481658830503[,[4.649481658830503;5.090963051959028[,[5.090963051959028;5.5324444450875525)>)
X4:Discretized(<(-3.3879472098791785;-2.9335217298098595[,[-2.9335217298098595;-2.479096249740541[,[-2.479096249740541;-2.024670769671222[,[-2.024670769671222;-1.570245289601903[,[-1.570245289601903;-1.1158198095325842[,[-1.1158198095325842;-0.6613943294632652[,[-0.6613943294632652;-0.20696884939394655[,[-0.20696884939394655;0.2474566306753725[,[0.2474566306753725;0.7018821107446911[,[0.7018821107446911;1.1563075908140101[,[1.1563075908140101;1.6107330708833292[,[1.6107330708833292;2.065158550952648[,[2.065158550952648;2.5195840310219673[,[2.5195840310219673;2.9740095110912854[,[2.9740095110912854;3.4284349911606045[,[3.4284349911606045;3.8828604712299235[,[3.8828604712299235;4.3372859512992425[,[4.3372859512992425;4.791711431368561[,[4.791711431368561;5.246136911437881[,[5.246136911437881;5.700562391507199)>)
Y:Discretized(<(-16.678566131823267;-13.748218943357312[,[-13.748218943357312;-10.817871754891357[,[-10.817871754891357;-7.8875245664254[,[-7.8875245664254;-4.957177377959445[,[-4.957177377959445;-2.0268301894934897[,[-2.0268301894934897;0.9035169989724672[,[0.9035169989724672;3.8338641874384223[,[3.8338641874384223;6.764211375904377[,[6.764211375904377;9.694558564370332[,[9.694558564370332;12.624905752836288[,[12.624905752836288;15.555252941302243[,[15.555252941302243;18.4856001297682[,[18.4856001297682;21.415947318234153[,[21.415947318234153;24.34629450670011[,[24.34629450670011;27.276641695166063[,[27.276641695166063;30.20698888363202[,[30.20698888363202;33.13733607209798[,[33.13733607209798;36.06768326056393[,[36.06768326056393;38.99803044902989[,[38.99803044902989;41.92837763749584)>)

We are now prepared to instantiate the CausalEffectEstimation object for both datasets.

linear_cee = gcee.CausalEffectEstimation(linear_data, disc_causal_model)
print(linear_cee)

<pyagrum.causalEffectEstimation._CausalEffectEstimation.CausalEffectEstimation object at 0x117750ad0>

 Dataframe  : <pandas.DataFrame object at 0x10fed7620>
  - shape    : (100000, 6)
  - columns  : Index(['T', 'X1', 'X2', 'X3', 'X4', 'Y'], dtype='str')
  - memory usage  : 4.800132 MB
 Causal Model  : <pyagrum.pyagrum.CausalModel; proxy of <Swig Object of type 'gum::CausalModel< double > *' at 0x1176e49b0> >
  - names    : {'X2', 'X3', 'X1', 'X4', 'Y', 'T'}
  - causal BN  : {0<T>,1<X1>,2<X2>,3<X3>,4<X4>,5<Y>} , {2->5,3->5,1->5,4->5,0->5}
  - observ. BN  : BN{nodes: 6, arcs: 5, domainSize: 10^6.80618, dim: 6080077, mem: 48Mo 848Ko 656o}

non_linear_cee = gcee.CausalEffectEstimation(non_linear_data, disc_causal_model)
print(non_linear_cee)

<pyagrum.causalEffectEstimation._CausalEffectEstimation.CausalEffectEstimation object at 0x1168c6850>

 Dataframe  : <pandas.DataFrame object at 0x1168c7110>
  - shape    : (100000, 6)
  - columns  : Index(['T', 'X1', 'X2', 'X3', 'X4', 'Y'], dtype='str')
  - memory usage  : 4.800132 MB
 Causal Model  : <pyagrum.pyagrum.CausalModel; proxy of <Swig Object of type 'gum::CausalModel< double > *' at 0x1176e49b0> >
  - names    : {'X2', 'X3', 'X1', 'X4', 'Y', 'T'}
  - causal BN  : {0<T>,1<X1>,2<X2>,3<X3>,4<X4>,5<Y>} , {2->5,3->5,1->5,4->5,0->5}
  - observ. BN  : BN{nodes: 6, arcs: 5, domainSize: 10^6.80618, dim: 6080077, mem: 48Mo 848Ko 656o}

1.3 Causal Indentification

The subsequent step involves identifying the causal criterion to be used for estimation. This is crucial, as most estimators rely on strong assumptions regarding the underlying causal structure of the data-generating process. Incorrect specification of the adjustment may compromise the guarantee of asymptotic normality.

linear_cee.identifyAdjustmentSet(intervention="T", outcome="Y")
non_linear_cee.identifyAdjustmentSet(intervention="T", outcome="Y", verbose=False)

Randomized Controlled Trial adjustment found.

Supported estimators include:
- CausalModelEstimator
- DM
If the outcome variable is a cause of other covariates in the causal graph,
Backdoor estimators may also be used.





'Randomized Controlled Trial'

Consistent with the data generation, the adjustment identified is the Randomized Control Trial adjustment. This yields a list of the various estimators supported by this adjustment.

1.4 Causal Effect Estimation

Once the adjustment is identified, we can proceed to estimation using the supported estimators. First, the estimator must be fitted to the data.

linear_cee.fitDM()
non_linear_cee.fitDM()

The estimation is obtained by calling the estimateCausalEffect method.

linear_tau_hat = linear_cee.estimateCausalEffect()
non_linear_tau_hat = non_linear_cee.estimateCausalEffect()

print(f"ACE linear = {linear_tau_hat}, MAPE = {abs((linear_tau_hat - 10) * 10)} %")
print(f"ACE non linear = {non_linear_tau_hat}, MAPE = {abs((non_linear_tau_hat - 10) * 10)} %")

ACE linear = 10.015848844809387, MAPE = 0.15848844809386975 %
ACE non linear = 10.02256525004448, MAPE = 0.22565250044479512 %

The difference-in-means estimator, which is the simplest estimator for the ACE, yields mostly accurate results. This is expected in an RCT environment where the treatment is independent of confounders, making intervention equivalent to observation.

1.5 User specified adjustment

If the user wish to apply an alternative adjustment, they may specify their own set of variables for each component of the adjustment. However, please note that such custom adjustments do not guarantee unbiased estimation and may not ensure an error-free estimation process.

linear_cee.useBackdoorAdjustment(intervention="T", outcome="Y", confounders={"X1", "X2", "X3", "X4"})
non_linear_cee.useBackdoorAdjustment(intervention="T", outcome="Y", confounders={"X1", "X2", "X3", "X4"})

linear_cee.fitSLearner()
non_linear_cee.fitSLearner()

linear_tau_hat = linear_cee.estimateCausalEffect()
non_linear_tau_hat = non_linear_cee.estimateCausalEffect()

print(f"ACE linear = {linear_tau_hat}, MAPE = {abs((linear_tau_hat - 10) * 10)} %")
print(f"ACE non linear = {non_linear_tau_hat}, MAPE = {abs((non_linear_tau_hat - 10) * 10)} %")

ACE linear = 10.038921959040907, MAPE = 0.3892195904090734 %
ACE non linear = 10.011560423132721, MAPE = 0.11560423132721098 %

In this case, RCT adjustment also supports Backdoor adjustment, we thus get mostly accurate estimations. Let’s see how the estimation would be if we specify an uncompatible adjustment.

1.6 Incorrectly specified adjustment

We will use the frontdoor adjustment to illustrate the behaviours of incorrect specification.

linear_cee.useFrontdoorAdjustment(intervention="T", outcome="Y", mediators={"X1", "X2", "X3", "X4"})
non_linear_cee.useFrontdoorAdjustment(intervention="T", outcome="Y", mediators={"X1", "X2", "X3", "X4"})

linear_cee.fitSimplePlugIn()
non_linear_cee.fitSimplePlugIn()

linear_tau_hat = linear_cee.estimateCausalEffect()
non_linear_tau_hat = non_linear_cee.estimateCausalEffect()

print(f"ACE linear = {linear_tau_hat}, MAPE = {abs((linear_tau_hat - 10) * 10)} %")
print(f"ACE non linear = {non_linear_tau_hat}, MAPE = {abs((non_linear_tau_hat - 10) * 10)} %")

ACE linear = -0.023219624375617063, MAPE = 100.23219624375616 %
ACE non linear = 0.010622858856612984, MAPE = 99.89377141143387 %

As anticipated, using an incorrect adjustment set results in a heavily biased estimation. In this case, the ACE is close to zero, indicating that the estimator incorrectly predicts no causal effect of $T$ on $Y$ .

1.7 Causal Bayesian Network Estimation

To fully utilize the causal model within the estimation process, we will now use the Conditional Probability Tables of the Causal Bayesian Network via the CausalBNEstimator, rather than relying solely on the underlying causal graph. The procedure will follow the same methodology as previously applied.

linear_cee.identifyAdjustmentSet(intervention="T", outcome="Y", verbose=False)
non_linear_cee.identifyAdjustmentSet(intervention="T", outcome="Y", verbose=False)

'Randomized Controlled Trial'

linear_cee.fitCausalBNEstimator()
non_linear_cee.fitCausalBNEstimator()

linear_tau_hat = linear_cee.estimateCausalEffect()
non_linear_tau_hat = non_linear_cee.estimateCausalEffect()

print(f"ACE linear = {linear_tau_hat}, MAPE = {abs((linear_tau_hat - 10) * 10)} %")
print(f"ACE non linear = {non_linear_tau_hat}, MAPE = {abs((non_linear_tau_hat - 10) * 10)} %")

ACE linear = 9.210563487180638, MAPE = 7.894365128193623 %
ACE non linear = 9.136945329244515, MAPE = 8.63054670755485 %

2. ACE estimations from real RCT data

In this section, we will estimate the ACE under real-world RCT conditions.

2.1 Dataset

The data used in this notbook come from the Tennessee Student/Teacher Achievement Ratio (STAR) trial. This randomized controlled trial was designed to assess the effects of smaller class sizes in primary schools (T) on students’ academic performance (Y).

The covariates in this study include:

gender
age
g1freelunch being the number of lunchs provided to the child per day
g1surban the localisation of the school (inner city or rural)
ethnicity

## Preprocessing

## Load data - read everything as a string and then cast
star_df = pd.read_csv("./res/STAR_data.csv", sep=",", dtype=str)
star_df = star_df.rename(columns={"race": "ethnicity"})

## Fill na
star_df = star_df.fillna({"g1freelunch": 0, "g1surban": 0})
drop_star_l = [
  "g1tlistss",
  "g1treadss",
  "g1tmathss",
  "g1classtype",
  "birthyear",
  "birthmonth",
  "birthday",
  "gender",
  "ethnicity",
  "g1freelunch",
  "g1surban",
]
star_df = star_df.dropna(subset=drop_star_l, how="any")

## Cast value types before processing
star_df["gender"] = star_df["gender"].astype(int)
star_df["ethnicity"] = star_df["ethnicity"].astype(int)

star_df["g1freelunch"] = star_df["g1freelunch"].astype(int)
star_df["g1surban"] = star_df["g1surban"].astype(int)
star_df["g1classtype"] = star_df["g1classtype"].astype(int)

## Keep only class type 1 and 2 (in the initial trial,
## 3 class types where attributed and the third one was big classes
## but with a teaching assistant)
star_df = star_df[~(star_df["g1classtype"] == 3)].reset_index(drop=True)

## Compute the outcome
star_df["Y"] = (
  star_df["g1tlistss"].astype(int) + star_df["g1treadss"].astype(int) + star_df["g1tmathss"].astype(int)
) / 3

## Compute the treatment
star_df["T"] = star_df["g1classtype"].apply(lambda x: 0 if x == 2 else 1)

## Transform date to obtain age (Notice: if na --> date is NaT)
star_df["date"] = pd.to_datetime(
  star_df["birthyear"] + "/" + star_df["birthmonth"] + "/" + star_df["birthday"], yearfirst=True, errors="coerce"
)
star_df["age"] = np.datetime64("1985-01-01") - star_df["date"]
star_df["age"] = star_df["age"].dt.days / 365.25

## Keep only covariates we consider predictive of the outcome
star_covariates_l = ["gender", "ethnicity", "age", "g1freelunch", "g1surban"]
star_df = star_df[["Y", "T"] + star_covariates_l]

## Map numerical to categorical
star_df["gender"] = star_df["gender"].apply(lambda x: "Girl" if x == 2 else "Boy").astype("category")
star_df["ethnicity"] = (
  star_df["ethnicity"]
  .map({1: "White", 2: "Black", 3: "Asian", 4: "Hispanic", 5: "Nat_American", 6: "Other"})
  .astype("category")
)
star_df["g1surban"] = (
  star_df["g1surban"].map({1: "Inner_city", 2: "Suburban", 3: "Rural", 4: "Urban"}).astype("category")
)

star_df.describe()

	Y	T	age	g1freelunch
count	4215.000000	4215.000000	4215.000000	4215.000000
mean	540.095848	0.428233	4.879872	1.471886
std	39.267221	0.494881	0.465104	0.534171
min	439.333333	0.000000	3.129363	0.000000
25%	511.333333	0.000000	4.525667	1.000000
50%	537.333333	0.000000	4.818617	1.000000
75%	566.000000	1.000000	5.111567	2.000000
max	670.666667	1.000000	7.225188	2.000000

It appears that there are more units in the control group. However, the control and treatment groups appear to be similar in distribution, indicating that the ignorability assumption is likely satisfied.

We will explore how the CausalEffectEstimation module can estimate the causal effect of $T$ on $Y$ in both of the given datasets.

2.2 Structure Learning and Setup

In the absence of a predefined causal structure, structure learning is utilized to uncover the underlying relationships between the variables in the dataset. To facilitate this process, a slice order will be imposed on the variables. This approach will serve as the foundation for deriving the necessary causal structure for subsequent analysis.

To enable the application of structure learning algorithms, the variables will first be discretized using the discretizer module. Following this, the causal structure will be derived using gum.BNLearner.

discretizer = disc.Discretizer(defaultDiscretizationMethod="uniform")
discretizer.setDiscretizationParameters("age", "uniform", 24)
discretizer.setDiscretizationParameters("Y", "uniform", 30)

template = discretizer.discretizedTemplate(star_df)

learner = gum.BNLearner(star_df, template)
learner.useNMLCorrection()
learner.useSmoothingPrior(1e-6)
learner.setSliceOrder(
  [
    ["T", "ethnicity", "gender", "age"],
    [
      "g1surban",
      "g1freelunch",
    ],
    ["Y"],
  ]
)
bn = learner.learnBN()

print(learner)

gnb.sideBySide(gexpl.getInformation(bn, size="50"), gnb.getInference(bn, size="50"))

Filename               : /var/folders/r1/pj4vdx_n4_d_xpsb04kzf97r0000gp/T/tmp_p1hx6gx.csv
Size                   : (4215,7)
Variables              : Y[30], T[2], gender[2], ethnicity[6], age[24], g1freelunch[3], g1surban[4]
Induced types          : False
Missing values         : False
Algorithm              : MIIC
Correction             : NML
Prior                  : Smoothing
Prior weight           : 0.000001
Constraint Slice Order : {ethnicity:0, T:0, g1surban:1, age:0, gender:0, g1freelunch:1, Y:2}

This initial approach appears promising, as the inferred causal relationships are somewhat consistent with what might be expected from an non-expert perspective.

Now given the causal structure, we are set to instanciate the CausalEffectEstimation class to perform estimation.

causal_model = gum.CausalModel(bn)
cee = gcee.CausalEffectEstimation(star_df, causal_model)

2.3 Causal Identification

The next step involves formal causal identification. As expected, we identify the RCT adjustment, consistent with the experimental design.

cee.identifyAdjustmentSet(intervention="T", outcome="Y")

Randomized Controlled Trial adjustment found.

Supported estimators include:
- CausalModelEstimator
- DM
If the outcome variable is a cause of other covariates in the causal graph,
Backdoor estimators may also be used.





'Randomized Controlled Trial'

2.4 Causal Effect Estimation

Once the ajustment identified, we can use the appropiate estimators for estimation.

cee.fitDM()
tau_hat = cee.estimateCausalEffect()

print(f"ACE = {tau_hat}")

ACE = 12.814738911047016

cee.fitCausalBNEstimator()
tau_hat = cee.estimateCausalEffect()

print(f"ACE = {tau_hat}")

ACE = 11.51523574877787

Let’s evaluate how the backdoor adjustment estimators compare to the previously obtained estimates. In this analysis, we control for the g1freelunch variable.

cee.useBackdoorAdjustment(intervention="T", outcome="Y", confounders={"g1freelunch"})

cee.fitSLearner()
tau_hat = cee.estimateCausalEffect()

print(f"ACE = {tau_hat}")

ACE = 11.616979201549668

cee.fitTLearner()
tau_hat = cee.estimateCausalEffect()

print(f"ACE = {tau_hat}")

ACE = 11.61670551692447

cee.fitIPW()
tau_hat = cee.estimateCausalEffect()

print(f"ACE = {tau_hat}")

ACE = 11.77382443916544

cee.fitPStratification()
tau_hat = cee.estimateCausalEffect()

print(f"ACE = {tau_hat}")

ACE = 10.906438153605961

3. ACE estimations from generated observational data

Next, we will estimate the ACE using generated data that does not adhere to RCT conditions. The generative model follows the specification from Lunceford & Davidian (2004)

3.1 Dataset

Consider the following generative model:

The outcome variable $Y$ is generated according the following equation: