From a7de9c135b0053c1166d2ddcebce25da16f4f112 Mon Sep 17 00:00:00 2001
From: sami bh <sami.boumaiza@ensta-paris.fr>
Date: Tue, 29 Aug 2023 17:25:20 +0200
Subject: [PATCH 1/2] added normalization

---
 src/benchmark_mediation.py | 54 +++++++++++++++++++++++++++-----------
 1 file changed, 38 insertions(+), 16 deletions(-)

diff --git a/src/benchmark_mediation.py b/src/benchmark_mediation.py
index 0c656f3..5e6eb31 100644
--- a/src/benchmark_mediation.py
+++ b/src/benchmark_mediation.py
@@ -628,6 +628,7 @@ def multiply_robust_efficient(
     forest=False,
     crossfit=0,
     trim=0.01,
+    normalized=True,
     regularization=True,
     calibration=True,
     calib_method="sigmoid",
@@ -668,6 +669,9 @@ def multiply_robust_efficient(
     trim : float, default=0.01
         Limit to trim p_x and f_mtx for numerical stability (min=trim, max=1-trim)
 
+    normalized : boolean, default=True
+        Normalizes the inverse probability-based weights.
+
     regularization : boolean, default=True
         Whether to use regularized models (logistic or linear regression).
         If True, cross-validation is used to chose among 8 potential
@@ -692,7 +696,7 @@ def multiply_robust_efficient(
         Direct effect on the exposed.
     direct0 : float
         Direct effect on the unexposed,
-    indirect1 : float 
+    indirect1 : float
         Indirect effect on the exposed.
     indirect0 : float
         Indirect effect on the unexposed.
@@ -898,27 +902,45 @@ def multiply_robust_efficient(
     f_m0x_trim = f_m0x != np.clip(f_m0x, trim, 1 - trim)
     f_m1x_trim = f_m1x != np.clip(f_m1x, trim, 1 - trim)
     trimmed = p_x_trim + f_m0x_trim + f_m1x_trim
-   
+
     var_name = ["t", "y", "p_x", "f_m0x", "f_m1x", "mu_t1", "mu_t0"]
     var_name += ["E_mu_t1_t1", "E_mu_t0_t0", "E_mu_t1_t0", "E_mu_t0_t1"]
-    
+
     for var in var_name:
         exec(f"{var} = {var}[~trimmed]")
     n_discarded += np.sum(trimmed)
 
-    # ytmt computing
-    y1m1 = t / p_x * (y - E_mu_t1_t1) + E_mu_t1_t1
-    y0m0 = (1 - t) / (1 - p_x) * (y - E_mu_t0_t0) + E_mu_t0_t0
-    y1m0 = (
-        (t / p_x) * (f_m0x / f_m1x) * (y - mu_t1)
-        + (1 - t) / (1 - p_x) * (mu_t1 - E_mu_t1_t0)
-        + E_mu_t1_t0
-    )
-    y0m1 = (
-        (1 - t) / (1 - p_x) * (f_m1x / f_m0x) * (y - mu_t0)
-        + t / p_x * (mu_t0 - E_mu_t0_t1)
-        + E_mu_t0_t1
-    )
+    # score computing
+    if normalized:
+        sumscore1 = np.mean(t / p_x)
+        sumscore2 = np.mean((1 - t) / (1 - p_x))
+        sumscore3 = np.mean((t / p_x) * (f_m0x / f_m1x))
+        sumscore4 = np.mean((1 - t) / (1 - p_x) * (f_m1x / f_m0x))
+        y1m1 = (t / p_x * (y - E_mu_t1_t1)) / sumscore1 + E_mu_t1_t1
+        y0m0 = ((1 - t) / (1 - p_x) * (y - E_mu_t0_t0)) / sumscore2 + E_mu_t0_t0
+        y1m0 = (
+            ((t / p_x) * (f_m0x / f_m1x) * (y - mu_t1)) / sumscore3
+            + ((1 - t) / (1 - p_x) * (mu_t1 - E_mu_t1_t0)) / sumscore2
+            + E_mu_t1_t0
+        )
+        y0m1 = (
+            ((1 - t) / (1 - p_x) * (f_m1x / f_m0x) * (y - mu_t0)) / sumscore4
+            + t / p_x * (mu_t0 - E_mu_t0_t1) / sumscore1
+            + E_mu_t0_t1
+        )
+    else:
+        y1m1 = t / p_x * (y - E_mu_t1_t1) + E_mu_t1_t1
+        y0m0 = (1 - t) / (1 - p_x) * (y - E_mu_t0_t0) + E_mu_t0_t0
+        y1m0 = (
+            (t / p_x) * (f_m0x / f_m1x) * (y - mu_t1)
+            + (1 - t) / (1 - p_x) * (mu_t1 - E_mu_t1_t0)
+            + E_mu_t1_t0
+        )
+        y0m1 = (
+            (1 - t) / (1 - p_x) * (f_m1x / f_m0x) * (y - mu_t0)
+            + t / p_x * (mu_t0 - E_mu_t0_t1)
+            + E_mu_t0_t1
+        )
 
     # effects computing
     total = np.mean(y1m1 - y0m0)

From 919ff9ad08b8c9f146b8d997998f7f3a22de1aec Mon Sep 17 00:00:00 2001
From: sami bh <sami.boumaiza@ensta-paris.fr>
Date: Thu, 31 Aug 2023 11:34:57 +0200
Subject: [PATCH 2/2] normalization documented

---
 src/benchmark_mediation.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/benchmark_mediation.py b/src/benchmark_mediation.py
index 5e6eb31..30d008d 100644
--- a/src/benchmark_mediation.py
+++ b/src/benchmark_mediation.py
@@ -670,7 +670,9 @@ def multiply_robust_efficient(
         Limit to trim p_x and f_mtx for numerical stability (min=trim, max=1-trim)
 
     normalized : boolean, default=True
-        Normalizes the inverse probability-based weights.
+        Normalizes the inverse probability-based weights so they add up to 1, as
+        described in "Identifying causal mechanisms (primarily) based on inverse probability weighting",
+        Huber (2014), https://doi.org/10.1002/jae.2341
 
     regularization : boolean, default=True
         Whether to use regularized models (logistic or linear regression).