Merge pull request #6 from YuminosukeSato/feat/summary-r-compat

feat: R-compatible summary format with full CausalImpactResults

Author: YuminosukeSato

python/causal_impact/analysis.py

@@ -11,21 +11,34 @@
 class CausalImpactResults:
     """Results of causal impact analysis."""
 
-    point_effects: np.ndarray  # (T_post,) mean effect per time point
-    point_effect_lower: np.ndarray  # (T_post,) lower CI per time point
-    point_effect_upper: np.ndarray  # (T_post,) upper CI per time point
-    ci_lower: float  # lower CI bound on average effect
-    ci_upper: float  # upper CI bound on average effect
-    point_effect_mean: float  # mean of point effects across time
-    cumulative_effect: np.ndarray  # (T_post,) cumulative point effects
-    cumulative_effect_lower: np.ndarray  # (T_post,) lower cumulative CI
-    cumulative_effect_upper: np.ndarray  # (T_post,) upper cumulative CI
-    cumulative_effect_total: float  # total cumulative effect
-    relative_effect_mean: float  # relative effect (effect / predicted)
-    p_value: float  # Bayesian one-sided tail probability
-    predictions_mean: np.ndarray  # (T_post,) mean counterfactual
-    predictions_lower: np.ndarray  # (T_post,) lower CI counterfactual
-    predictions_upper: np.ndarray  # (T_post,) upper CI counterfactual
+    actual: np.ndarray
+    point_effects: np.ndarray
+    point_effect_lower: np.ndarray
+    point_effect_upper: np.ndarray
+    ci_lower: float
+    ci_upper: float
+    point_effect_mean: float
+    average_effect_sd: float
+    cumulative_effect: np.ndarray
+    cumulative_effect_lower: np.ndarray
+    cumulative_effect_upper: np.ndarray
+    cumulative_effect_total: float
+    cumulative_effect_sd: float
+    relative_effect_mean: float
+    relative_effect_sd: float
+    relative_effect_lower: float
+    relative_effect_upper: float
+    p_value: float
+    predictions_mean: np.ndarray
+    predictions_sd: np.ndarray
+    predictions_lower: np.ndarray
+    predictions_upper: np.ndarray
+    average_prediction_sd: float
+    average_prediction_lower: float
+    average_prediction_upper: float
+    cumulative_prediction_sd: float
+    cumulative_prediction_lower: float
+    cumulative_prediction_upper: float
 
 
 class CausalAnalysis:
@@ -42,28 +55,19 @@ def compute_effects(
 
         n_samples = predictions.shape[0]
 
-        # Effect per sample per time point: observed - counterfactual
-        # predictions shape: (n_samples, t_post)
-        effects = y_post[np.newaxis, :] - predictions  # (n_samples, t_post)
+        effects = y_post[np.newaxis, :] - predictions
+        avg_effects = effects.mean(axis=1)
+        point_effects = effects.mean(axis=0)
 
-        # Average effect across time for each sample
-        avg_effects = effects.mean(axis=1)  # (n_samples,)
-
-        # Point effects: mean across samples at each time point
-        point_effects = effects.mean(axis=0)  # (t_post,)
-
-        # Summary-table CI on average effect uses sample-average quantiles.
         lower_q = alpha / 2
         upper_q = 1 - alpha / 2
         point_effect_lower = np.percentile(effects, 100 * lower_q, axis=0)
         point_effect_upper = np.percentile(effects, 100 * upper_q, axis=0)
         ci_lower = float(np.percentile(avg_effects, 100 * lower_q))
         ci_upper = float(np.percentile(avg_effects, 100 * upper_q))
 
-        # Mean effect
         point_effect_mean = float(avg_effects.mean())
 
-        # Cumulative effect
         cumulative_effect = np.cumsum(point_effects)
         cum_effects_samples = np.cumsum(effects, axis=1)
         cumulative_effect_lower = np.percentile(
@@ -78,40 +82,107 @@ def compute_effects(
         )
         cumulative_effect_total = float(cumulative_effect[-1])
 
-        # Relative effect
-        pred_mean_total = predictions.mean()
-        if abs(pred_mean_total) > 1e-10:
-            relative_effect_mean = point_effect_mean / pred_mean_total
+        actual = y_post.copy()
+
+        if n_samples == 1:
+            predictions_sd_arr = np.zeros(predictions.shape[1])
+        else:
+            predictions_sd_arr = np.std(predictions, axis=0, ddof=1)
+
+        avg_pred_per_sample = predictions.mean(axis=1)
+        cum_pred_per_sample = predictions.sum(axis=1)
+
+        if n_samples == 1:
+            average_prediction_sd = 0.0
+            cumulative_prediction_sd = 0.0
         else:
-            relative_effect_mean = 0.0
+            average_prediction_sd = float(np.std(avg_pred_per_sample, ddof=1))
+            cumulative_prediction_sd = float(np.std(cum_pred_per_sample, ddof=1))
+
+        average_prediction_lower = float(
+            np.percentile(avg_pred_per_sample, 100 * lower_q)
+        )
+        average_prediction_upper = float(
+            np.percentile(avg_pred_per_sample, 100 * upper_q)
+        )
+        cumulative_prediction_lower = float(
+            np.percentile(cum_pred_per_sample, 100 * lower_q)
+        )
+        cumulative_prediction_upper = float(
+            np.percentile(cum_pred_per_sample, 100 * upper_q)
+        )
+
+        cum_effects_per_sample = effects.sum(axis=1)
+
+        if n_samples == 1:
+            average_effect_sd = 0.0
+            cumulative_effect_sd = 0.0
+        else:
+            average_effect_sd = float(np.std(avg_effects, ddof=1))
+            cumulative_effect_sd = float(np.std(cum_effects_per_sample, ddof=1))
+
+        avg_pred_per_sample_safe = np.where(
+            np.abs(avg_pred_per_sample) > 1e-10,
+            avg_pred_per_sample,
+            np.nan,
+        )
+        rel_effects_per_sample = np.where(
+            np.abs(avg_pred_per_sample) > 1e-10,
+            avg_effects / avg_pred_per_sample_safe,
+            0.0,
+        )
+
+        relative_effect_mean = float(rel_effects_per_sample.mean())
+
+        if n_samples == 1:
+            relative_effect_sd = 0.0
+        else:
+            relative_effect_sd = float(np.std(rel_effects_per_sample, ddof=1))
+
+        relative_effect_lower = float(
+            np.percentile(rel_effects_per_sample, 100 * lower_q)
+        )
+        relative_effect_upper = float(
+            np.percentile(rel_effects_per_sample, 100 * upper_q)
+        )
 
-        # p-value: proportion of samples where average effect has opposite sign
         if point_effect_mean >= 0:
             p_value = float(np.mean(avg_effects < 0))
         else:
             p_value = float(np.mean(avg_effects > 0))
-        # Ensure minimum p-value of 1/n_samples
         p_value = max(p_value, 1.0 / n_samples)
 
-        # Counterfactual prediction summaries
         predictions_mean = predictions.mean(axis=0)
         predictions_lower = np.percentile(predictions, 100 * lower_q, axis=0)
         predictions_upper = np.percentile(predictions, 100 * upper_q, axis=0)
 
         return CausalImpactResults(
+            actual=actual,
             point_effects=point_effects,
             point_effect_lower=point_effect_lower,
             point_effect_upper=point_effect_upper,
             ci_lower=ci_lower,
             ci_upper=ci_upper,
             point_effect_mean=point_effect_mean,
+            average_effect_sd=average_effect_sd,
             cumulative_effect=cumulative_effect,
             cumulative_effect_lower=cumulative_effect_lower,
             cumulative_effect_upper=cumulative_effect_upper,
             cumulative_effect_total=cumulative_effect_total,
+            cumulative_effect_sd=cumulative_effect_sd,
             relative_effect_mean=relative_effect_mean,
+            relative_effect_sd=relative_effect_sd,
+            relative_effect_lower=relative_effect_lower,
+            relative_effect_upper=relative_effect_upper,
             p_value=p_value,
             predictions_mean=predictions_mean,
+            predictions_sd=predictions_sd_arr,
             predictions_lower=predictions_lower,
             predictions_upper=predictions_upper,
+            average_prediction_sd=average_prediction_sd,
+            average_prediction_lower=average_prediction_lower,
+            average_prediction_upper=average_prediction_upper,
+            cumulative_prediction_sd=cumulative_prediction_sd,
+            cumulative_prediction_lower=cumulative_prediction_lower,
+            cumulative_prediction_upper=cumulative_prediction_upper,
         )

python/causal_impact/summary.py

@@ -6,42 +6,83 @@
 
 
 class SummaryFormatter:
-    """Format CausalImpact results as text summary or natural language report."""
+    """Format CausalImpact results as text summary or report."""
 
     @staticmethod
     def summary(results: CausalImpactResults, digits: int = 2) -> str:
         fmt = f".{digits}f"
 
-        avg_effect = format(results.point_effect_mean, fmt)
-        avg_ci = f"[{format(results.ci_lower, fmt)}, {format(results.ci_upper, fmt)}]"
-        cum_effect = format(results.cumulative_effect_total, fmt)
-        cum_ci = (
+        avg_actual = format(results.actual.mean(), fmt)
+        cum_actual = format(results.actual.sum(), fmt)
+
+        avg_pred = format(results.predictions_mean.mean(), fmt)
+        avg_pred_sd = format(results.average_prediction_sd, fmt)
+        cum_pred = format(results.predictions_mean.sum(), fmt)
+        cum_pred_sd = format(results.cumulative_prediction_sd, fmt)
+
+        avg_pred_ci = (
+            f"[{format(results.average_prediction_lower, fmt)}, "
+            f"{format(results.average_prediction_upper, fmt)}]"
+        )
+        cum_pred_ci = (
+            f"[{format(results.cumulative_prediction_lower, fmt)}, "
+            f"{format(results.cumulative_prediction_upper, fmt)}]"
+        )
+
+        avg_eff = format(results.point_effect_mean, fmt)
+        avg_eff_sd = format(results.average_effect_sd, fmt)
+        cum_eff = format(results.cumulative_effect_total, fmt)
+        cum_eff_sd = format(results.cumulative_effect_sd, fmt)
+
+        avg_eff_ci = (
+            f"[{format(results.ci_lower, fmt)}, {format(results.ci_upper, fmt)}]"
+        )
+        cum_eff_ci = (
             f"[{format(results.cumulative_effect_lower[-1], fmt)}, "
             f"{format(results.cumulative_effect_upper[-1], fmt)}]"
         )
-        rel_effect = format(results.relative_effect_mean * 100, fmt)
+
+        rel_m = format(results.relative_effect_mean * 100, fmt)
+        rel_sd = format(results.relative_effect_sd * 100, fmt)
+        rel_lo = format(results.relative_effect_lower * 100, fmt)
+        rel_hi = format(results.relative_effect_upper * 100, fmt)
+
         p_val = format(results.p_value, f".{max(digits, 3)}f")
+        prob = format((1 - results.p_value) * 100, fmt)
+
+        pred_row = (
+            f"Prediction (s.d.)        "
+            f"{avg_pred} ({avg_pred_sd})   "
+            f"{cum_pred} ({cum_pred_sd})"
+        )
+        eff_row = (
+            f"Absolute effect (s.d.)   "
+            f"{avg_eff} ({avg_eff_sd})    "
+            f"{cum_eff} ({cum_eff_sd})"
+        )
+        rel_row = f"Relative effect (s.d.)   {rel_m}% ({rel_sd}%) {rel_m}% ({rel_sd}%)"
+        rel_ci_row = (
+            f"95% CI                   [{rel_lo}%, {rel_hi}%] [{rel_lo}%, {rel_hi}%]"
+        )
 
         lines = [
             "Posterior inference {CausalImpact}",
             "",
             "                         Average        Cumulative",
-            "Actual                   -              -",
-            "Prediction (s.d.)        -              -",
-            f"95% CI                   {avg_ci}       {cum_ci}",
+            f"Actual                   {avg_actual}          {cum_actual}",
+            pred_row,
+            f"95% CI                   {avg_pred_ci}  {cum_pred_ci}",
+            "",
+            eff_row,
+            f"95% CI                   {avg_eff_ci}   {cum_eff_ci}",
             "",
-            f"Absolute effect (mean)   {avg_effect}           {cum_effect}",
-            f"Relative effect          {rel_effect}%",
+            rel_row,
+            rel_ci_row,
             "",
             f"Posterior tail-area probability p: {p_val}",
+            f"Posterior prob. of a causal effect: {prob}%",
         ]
 
-        if results.p_value < 0.05:
-            lines.append("Posterior prob. of a causal effect: "
-                         f"{format((1 - results.p_value) * 100, fmt)}%")
-        else:
-            lines.append("The effect is not statistically significant.")
-
         return "\n".join(lines)
 
     @staticmethod