Statistical Evidence

This section formalises the causal argument with statistical tests. The analysis uses vendor teams as a natural control group for the skill-based routing change — vendor teams were subject to the same self-serve automation and the same case mix shift, but were not subject to skill-based routing. Any difference in SLA trajectory between internal and vendor teams after Week 14 is therefore attributable to routing.

Test design: - Chi-square test — confirms internal and vendor SLA distributions were statistically equivalent at baseline and significantly different post-routing - Difference-in-Differences (DiD) — estimates the causal effect of routing by comparing internal vs vendor SLA change from baseline to both-changes period

Assumption: Parallel trends — in the absence of routing, internal and vendor SLA would have followed similar trajectories. This is supported by the baseline period where both groups tracked within 1–2pp of each other, and by the self-serve only period where both groups declined similarly.

Limitation: VendorA showed an independent performance dip around Week 14. This may slightly understate the true routing effect since the vendor control group is imperfect.

Setup and Load Data

Code 
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from statsmodels.stats.proportion import proportions_ztest
from scipy.stats import chi2_contingency
import warnings

warnings.filterwarnings("ignore")

cases = pd.read_csv("ops_cases.csv", parse_dates=["created_at", "resolved_at"])
period_order = ["Baseline", "Self-serve Only", "Both Changes"]
cases["period"] = pd.Categorical(cases["period"], categories=period_order, ordered=True)
cases.shape
(12841, 17)

1. Chi-Square Test - Baseline vs Both Changes

Setup

We test whether SLA outcome is independent of handler type (Internal/Vendor) within two distinct periods:

Period A: Baseline (Weeks 1–10) — before any operational changes Period B: Both Changes (Weeks 14–26) — after both self-serve and routing are live

The test is run separately for each period.

Period A — Baseline Period

\(H_{0}\): The proportion of cases resolved within SLA is the same for Internal and Vendor handler types during the baseline period.

P(SLA pass | Internal, Baseline) = P(SLA pass | Vendor, Baseline)

\(H_{1}\): The proportions differ.

P(SLA pass | Internal, Baseline) ≠ P(SLA pass | Vendor, Baseline)

Significance level: α = 0.05

Result: χ²(1) = 2.101, p = 0.147

Interpretation: We fail to reject H₀. There is no statistically significant difference in SLA pass rates between Internal and Vendor teams during the baseline period. This is a critical validation — it confirms the two groups were operationally comparable before any intervention, supporting the use of vendor teams as a control group in the DiD analysis that follows.

Period B — Both Changes

\(H_{0}\): The proportion of cases resolved within SLA is the same for Internal and Vendor handler types during the both-changes period.

P(SLA pass | Internal, Both Changes) = P(SLA pass | Vendor, Both Changes)

\(H_{1}\): The proportions differ.

P(SLA pass | Internal, Both Changes) ≠ P(SLA pass | Vendor, Both Changes)

Significance level: α = 0.05

Result: χ²(1) = 81.556, p < 0.001

Interpretation: We reject H₀. There is a significant difference in SLA rates between Internal and Vendor teams after both changes are live.

Groups that were statistically equivalent at baseline have diverged sharply. Since the only change that affected Internal teams specifically was skill-based routing, this divergence is attributable to the routing change.

Assumptions and Validity

  1. Independence of observations - each case is a unique event.
  2. Expected cell counts ≥ 5 - satisfied given sample sizes of thousands per cell.
  3. Mutual exclusivity - each case belongs to exactly one handler type and one SLA outcome.
Code 
for period in ["Baseline", "Both Changes"]:
    subset = cases[cases["period"] == period]
    ct = pd.crosstab(subset["initial_handler_type"], subset["resolved_within_sla"])
    chi2, p, dof, expected = chi2_contingency(ct)
    print(f"\n{'='*40}")
    print(f"Period: {period}")
    print(ct)
    sla_rates = ct.div(ct.sum(axis=1), axis=0)[1].rename("SLA rate")
    print(sla_rates.round(3))
    print(f"Chi2: {chi2:.3f}, p-value: {p:.4f}")
    if p > 0.05:
        print("No significant difference between groups (validates control group)")
    else:
        print("Significant difference detected")

========================================
Period: Baseline
resolved_within_sla     0     1
initial_handler_type           
Internal              553  2848
Vendor                417  1932
initial_handler_type
Internal    0.837
Vendor      0.822
Name: SLA rate, dtype: float64
Chi2: 2.101, p-value: 0.1472
No significant difference between groups (validates control group)

========================================
Period: Both Changes
resolved_within_sla      0     1
initial_handler_type            
Internal              1196  2290
Vendor                 523  1740
initial_handler_type
Internal    0.657
Vendor      0.769
Name: SLA rate, dtype: float64
Chi2: 81.556, p-value: 0.0000
Significant difference detected

2. Difference-in-Differences Estimate

Straightforward Approach

Code 
internal_baseline = cases[
    (cases["period"] == "Baseline") & (cases["initial_handler_type"] == "Internal")
]["resolved_within_sla"].mean()
internal_post = cases[
    (cases["period"] == "Both Changes") & (cases["initial_handler_type"] == "Internal")
]["resolved_within_sla"].mean()
vendor_baseline = cases[
    (cases["period"] == "Baseline") & (cases["initial_handler_type"] == "Vendor")
]["resolved_within_sla"].mean()
vendor_post = cases[
    (cases["period"] == "Both Changes") & (cases["initial_handler_type"] == "Vendor")
]["resolved_within_sla"].mean()

internal_change = internal_post - internal_baseline
vendor_change = vendor_post - vendor_baseline
did_estimate = internal_change - vendor_change

summary = pd.DataFrame(
    {
        "Baseline Period SLA": [f"{internal_baseline:.1%}", f"{vendor_baseline:.1%}"],
        "Both Changes Period SLA": [f"{internal_post:.1%}", f"{vendor_post:.1%}"],
        "Change": [f"{internal_change:+.1%}", f"{vendor_change:+.1%}"],
    },
    index=["Internal", "Vendor"],
)

print(summary)
print(f"\nDiD estimate (routing effect): {did_estimate:+.1%}")
         Baseline Period SLA Both Changes Period SLA  Change
Internal               83.7%                   65.7%  -18.0%
Vendor                 82.2%                   76.9%   -5.4%

DiD estimate (routing effect): -12.7%
Code 
# Contribution of Internal case SLA decline
internal_prop = (
    cases[cases["initial_handler_type"] == "Internal"].shape[0] / cases.shape[0]
)

print(f"Contribution of Internal Case SLA decline: {did_estimate*internal_prop:+.2%}")
Contribution of Internal Case SLA decline: -7.61%

Interpretation: After accounting for the -5.4% decline experienced by vendor teams (background + mix shift effect), skill-based routing is estimated to have caused an additional 12.7% decline in internal team SLA.

DiD Regression Analysis

Regression Equation for this problem statement:

\(SLA = β0 + β1·Treated + β2·Post + β3·(Treated × Post) + ε\)

  • β0 = baseline SLA for vendor (control) in pre-period
  • β1 = pre-existing gap between internal and vendor (should be near zero)
  • β2 = time trend affecting everyone (vendor’s change from pre to post)
  • β3 = the DiD estimate — the additional effect on internal teams beyond the background trend due to the treatment
Code 
# Restrict to Baseline and Both Changes only
did_data = cases[cases["period"].isin(["Baseline", "Both Changes"])].copy()

# Binary indicators for eaach group
did_data["treated"] = (did_data["initial_handler_type"] == "Internal").astype(int)
did_data["post"] = (did_data["period"] == "Both Changes").astype(int)
did_data["treated_post"] = did_data["treated"] * did_data["post"]

# Basic DiD — Linear Probability Model (LPM)
# LPM is standard for DiD with binary outcomes — coefficients are directly interpretable as pp changes
model1 = smf.ols("resolved_within_sla ~ treated + post + treated_post", data=did_data)
result1 = model1.fit(cov_type="HC3")  # HC3 = heteroscedasticity-robust standard errors

print("=== Model 1: Basic DiD ===")
print(result1.summary2().tables[1].round(4))
print(
    f"\nDiD estimate:      {result1.params['treated_post']:.4f} ({result1.params['treated_post']*100:.2f}pp)"
)
print(
    f"95% CI:            [{result1.conf_int().loc['treated_post', 0]*100:.2f}pp, "
    f"{result1.conf_int().loc['treated_post', 1]*100:.2f}pp]"
)
print(f"p-value:           {result1.pvalues['treated_post']}")
=== Model 1: Basic DiD ===
               Coef.  Std.Err.         z  P>|z|  [0.025  0.975]
Intercept     0.8225    0.0079  104.2779   0.00  0.8070  0.8379
treated       0.0149    0.0101    1.4757   0.14 -0.0049  0.0347
post         -0.0536    0.0119   -4.5160   0.00 -0.0768 -0.0303
treated_post -0.1269    0.0157   -8.0983   0.00 -0.1576 -0.0962

DiD estimate:      -0.1269 (-12.69pp)
95% CI:            [-15.76pp, -9.62pp]
p-value:           5.574586632611151e-16
Code 
# Add queue_name as a covariate
model2 = smf.ols(
    "resolved_within_sla ~ treated + post + treated_post + C(queue_name)", data=did_data
)
result2 = model2.fit(cov_type="HC3")

print("\n=== Model 2: DiD with Queue Controls ===")
print(result2.summary2().tables[1].round(4))
print(
    f"\nDiD estimate (controlled): {result2.params['treated_post']:.4f} "
    f"({result2.params['treated_post']*100:.2f}pp)"
)
print(
    f"95% CI:                    [{result2.conf_int().loc['treated_post', 0]*100:.2f}pp, "
    f"{result2.conf_int().loc['treated_post', 1]*100:.2f}pp]"
)

=== Model 2: DiD with Queue Controls ===
                                     Coef.  Std.Err.        z   P>|z|  [0.025  \
Intercept                           0.8327    0.0089  93.5151  0.0000  0.8153   
C(queue_name)[T.Complex_Cases]     -0.0304    0.0108  -2.8218  0.0048 -0.0515   
C(queue_name)[T.Technical_Support] -0.0102    0.0091  -1.1244  0.2608 -0.0280   
treated                             0.0145    0.0101   1.4420  0.1493 -0.0052   
post                               -0.0483    0.0120  -4.0324  0.0001 -0.0718   
treated_post                       -0.1266    0.0157  -8.0822  0.0000 -0.1573   

                                    0.975]  
Intercept                           0.8502  
C(queue_name)[T.Complex_Cases]     -0.0093  
C(queue_name)[T.Technical_Support]  0.0076  
treated                             0.0343  
post                               -0.0248  
treated_post                       -0.0959  

DiD estimate (controlled): -0.1266 (-12.66pp)
95% CI:                    [-15.73pp, -9.59pp]

Routing Changes Test Summary

  • Baseline chi-square (p = 0.147): Internal and vendor SLA were statistically equivalent before any intervention. This validates the parallel trends assumption required for a valid DiD estimate.

  • Both Changes chi-square (p < 0.001): After both changes are live, the two groups are statistically significantly different. Groups that were identical at baseline have diverged — and the only change that affected internal teams specifically was skill-based routing.

  • DiD estimate: −12.7pp. Internal SLA declined 18.0pp from baseline to the both-changes period. Vendor SLA declined 5.4pp over the same window, capturing the background effect (mix shift + general difficulty of remaining caseload). Stripping that out, skill-based routing is estimated to have caused an additional −12.7 percentage point decline in internal SLA.

  • Limitation: The vendor control group is imperfect. VendorA showed an independent performance dip around Week 14 unrelated to routing, which may slightly understate the true routing effect. The true causal estimate may be larger than −12.7pp.

Self- Serve Changes Assessment

“How much of the SLA decline after Week 11 can be explained purely by the change in case complexity mix caused by self-serve deflection?”

Code 
# Isolate the data for the self serve period - not considering the routing changes period
cases["period_selfserve"] = cases["week_number"].apply(
    lambda x: "pre" if x < 11 else "post"
)
ss_data = cases[cases["week_number"] < 14]

Proportions z test

Did the proportion of each queue change significantly after the self-serve period

\(H_{0}\): There is no difference between the two proportions

\(H_{1}\): There is significant difference between the two proportions

Code 
# Run the prortions z test
total_pre = len(ss_data[ss_data["period_selfserve"] == "pre"])
total_post = len(ss_data[ss_data["period_selfserve"] == "post"])

results = []
for queue in cases["queue_name"].unique():
    count_pre = len(
        ss_data[
            (ss_data["period_selfserve"] == "pre") & (ss_data["queue_name"] == queue)
        ]
    )
    count_post = len(
        ss_data[
            (ss_data["period_selfserve"] == "post") & (ss_data["queue_name"] == queue)
        ]
    )

    # Test: did this queue's share of total volume change?
    stat, p = proportions_ztest(
        count=[count_pre, count_post], nobs=[total_pre, total_post]
    )

    results.append(
        {
            "queue": queue,
            "share_pre": round(count_pre / total_pre, 3),
            "share_post": round(count_post / total_post, 3),
            "share_delta": round(
                (count_post / total_post) - (count_pre / total_pre), 3
            ),
            "z_stat": round(stat, 3),
            "p_value": round(p, 4),
            "significant": p < 0.05,
        }
    )

pd.DataFrame(results).sort_values("p_value")
queue share_pre share_post share_delta z_stat p_value significant
0 Basic_Support 0.405 0.189 -0.216 14.747 0.0000 True
2 Complex_Cases 0.198 0.369 0.171 -13.392 0.0000 True
1 Technical_Support 0.397 0.442 0.045 -2.991 0.0028 True

Quantify mechanical SLA impact from mix shift

If agent performance within each complexity bucket stayed constant, how much SLA decline would occur purely because the mix changed? Using counterfactuals

Code 
sla_table = (
    cases[cases["period"] == "Baseline"]
    .groupby(["queue_name"], as_index=False)
    .agg(sla_rate=("resolved_within_sla", "mean"))
    .round(3)
)

prop_baseline = (
    pd.DataFrame(
        cases[cases["period"] == "Baseline"]["queue_name"].value_counts(normalize=True)
    )
    .reset_index()
    .round(3)
)
prop_ss = (
    pd.DataFrame(
        cases[cases["period"] == "Self-serve Only"]["queue_name"].value_counts(
            normalize=True
        )
    )
    .reset_index()
    .round(3)
)
sla_table = pd.merge(sla_table, prop_ss, on="queue_name")

sla_table
queue_name sla_rate proportion
0 Basic_Support 0.890 0.199
1 Complex_Cases 0.753 0.360
2 Technical_Support 0.810 0.441
Code 
# Counterfactual SLA
print(
    "Expected SLA if only the workload mix changed.: ",
    (np.dot(sla_table["sla_rate"], sla_table["proportion"]).round(3)),
)
Expected SLA if only the workload mix changed.:  0.805

Interpret the Results

  • Actual Baseline SLA: 83.1%
  • Expected SLA from mix shift: 80.5%
  • Actual observed Weeks 11–13 SLA: 78.8%

Results

  • 2.6pp decline explained by complexity mix
  • Additional 1.7pp decline likely effects of solving harder cases and time spent increasing
Code 
# Run the prortions z test to see if the SLA decline is significant or not

# Split periods
pre = cases[cases["period"] == "Baseline"]
post = cases[cases["period"] == "Self-serve Only"]

# Number of SLA successes
sla_pre = pre["resolved_within_sla"].sum()
sla_post = post["resolved_within_sla"].sum()

# Total observations
n_pre = len(pre)
n_post = len(post)

# Two proportion z-test
stat, p = proportions_ztest(count=[sla_pre, sla_post], nobs=[n_pre, n_post])

print("Pre SLA:", round(sla_pre / n_pre, 3))
print("Post SLA:", round(sla_post / n_post, 3))
print("Delta:", round((sla_post / n_post) - (sla_pre / n_pre), 3))

print("Z-stat:", round(stat, 3))
print("P-value:", round(p, 5))
print("Significant:", p < 0.05)
Pre SLA: 0.831
Post SLA: 0.788
Delta: -0.044
Z-stat: 3.775
P-value: 0.00016
Significant: True

The SLA decline following Week 11 is statistically significant and unlikely due to random operational variation.

Sent Back Rate Analysis

Code 
reviews = pd.read_csv(
    "ops_reviews.csv", parse_dates=["intended_start_at", "actual_start_at", "end_at"]
)
reviews.shape
(17135, 11)
Code 
# Sent-back rate by skill bucket — requires routing config data
# As proxy, use queue_name × agent_id from ops_reviews to find which agent + case type pairings generate the most send-backs

sent_back_by_agent = (
    reviews[(reviews["outcome"] == "Sent Back")]
    .groupby(["agent_id", "queue_name", "handler_type"])
    .agg(sent_back_count=("review_id", "count"))
    .reset_index()
)

total_by_agent = (
    reviews.groupby(["agent_id", "queue_name", "handler_type"])
    .agg(total_reviews=("review_id", "count"))
    .reset_index()
)

agent_sent_back = sent_back_by_agent.merge(
    total_by_agent, on=["agent_id", "queue_name", "handler_type"]
)
agent_sent_back["sent_back_rate"] = (
    agent_sent_back["sent_back_count"] / agent_sent_back["total_reviews"]
)
Code 
agent_sent_back.sort_values(["sent_back_rate"], ascending=False)
agent_id queue_name handler_type sent_back_count total_reviews sent_back_rate
597 VC028 Complex_Cases Vendor 4 6 0.666667
97 AGT033 Complex_Cases Internal 23 39 0.589744
640 VC047 Complex_Cases Vendor 4 7 0.571429
172 AGT058 Complex_Cases Internal 22 39 0.564103
595 VC027 Complex_Cases Vendor 5 9 0.555556
... ... ... ... ... ... ...
373 VA031 Technical_Support Vendor 1 20 0.050000
535 VB050 Complex_Cases Vendor 1 21 0.047619
581 VC021 Technical_Support Vendor 1 23 0.043478
87 AGT030 Basic_Support Internal 1 23 0.043478
430 VB004 Technical_Support Vendor 1 24 0.041667

650 rows × 6 columns

Code 
# Merge period and handler type into reviews
reviews = reviews.merge(
    cases[["case_id", "period", "initial_handler_type"]], on="case_id", how="left"
)

# Per-review timing
reviews["handle_hours"] = (
    reviews["end_at"] - reviews["actual_start_at"]
).dt.total_seconds() / 3600

# Sort and shift to get time to next review pick-up within each case
reviews = reviews.sort_values(["case_id", "review_order"])
reviews["next_actual_start"] = reviews.groupby("case_id")["actual_start_at"].shift(-1)

reviews["requeue_wait_hours"] = (
    reviews["next_actual_start"] - reviews["end_at"]
).dt.total_seconds() / 3600

# Isolate sent-back reviews
sb = reviews[reviews["outcome"] == "Sent Back"].copy()
sb["total_cost_hours"] = sb["handle_hours"] + sb["requeue_wait_hours"]

# Isolate escalations
esc = reviews[reviews["outcome"] == "Escalated"].copy()
esc["total_cost_hours"] = esc["handle_hours"] + esc["requeue_wait_hours"]

sb.shape, esc.shape
((4082, 17), (212, 17))
Code 
# Time cost per sent-back event
cost_summary = (
    sb.groupby(["period", "handler_type"])
    .agg(
        avg_failed_review_hrs=("handle_hours", "mean"),
        avg_requeue_wait_hrs=("requeue_wait_hours", "mean"),
        avg_total_cost_hrs=("total_cost_hours", "mean"),
        num_sent_back_events=("review_id", "count"),
    )
    .round(2)
)

print("=== Time cost per send-back event ===")
cost_summary
=== Time cost per send-back event ===
avg_failed_review_hrs avg_requeue_wait_hrs avg_total_cost_hrs num_sent_back_events
period handler_type
Baseline Internal 0.71 29.81 30.52 565
Vendor 0.66 31.69 32.34 379
Self-serve Only Internal 0.86 33.19 34.05 183
Vendor 0.75 34.19 34.94 122
Both Changes Internal 0.78 28.78 29.56 2358
Vendor 0.77 36.03 36.80 475

1700 events can be prevented potentially (2358 - 650)

Code 
# Count send-back events per case
sb_counts = (
    reviews[reviews["outcome"] == "Sent Back"]
    .groupby("case_id")
    .size()
    .reset_index(name="num_sent_backs")
)

cases_reg = cases.merge(sb_counts, on="case_id", how="left")
cases_reg["num_sent_backs"] = cases_reg["num_sent_backs"].fillna(0)

# Run separately for Baseline and Both Changes
for period in ["Baseline", "Self-serve Only", "Both Changes"]:
    subset = cases_reg[cases_reg["period"] == period]

    model = smf.ols(
        "resolution_hours ~ num_sent_backs + C(queue_name) + C(initial_handler_type)",
        data=subset,
    ).fit(cov_type="HC3")

    coef = model.params["num_sent_backs"]
    pval = model.pvalues["num_sent_backs"]
    ci = model.conf_int().loc["num_sent_backs"]

    print(f"\n{period}")
    print(f"  Hours added per send-back event:  {coef:.2f}h")
    print(f"  95% CI:                           [{ci[0]:.2f}h, {ci[1]:.2f}h]")
    print(f"  p-value:                          {pval:.4f}")

Baseline
  Hours added per send-back event:  -0.28h
  95% CI:                           [-1.46h, 0.91h]
  p-value:                          0.6455

Self-serve Only
  Hours added per send-back event:  1.25h
  95% CI:                           [-1.22h, 3.72h]
  p-value:                          0.3222

Both Changes
  Hours added per send-back event:  3.95h
  95% CI:                           [2.97h, 4.92h]
  p-value:                          0.0000
Code 
# Count send-back events per case
esc_counts = (
    reviews[reviews["outcome"] == "Escalated"]
    .groupby("case_id")
    .size()
    .reset_index(name="num_escalations")
)

cases_reg = cases.merge(esc_counts, on="case_id", how="left")
cases_reg["num_escalations"] = cases_reg["num_escalations"].fillna(0)

# Run separately for Baseline and Both Changes
for period in ["Baseline", "Both Changes"]:
    subset = cases_reg[cases_reg["period"] == period]

    model = smf.ols(
        "resolution_hours ~ num_escalations + C(queue_name) + C(initial_handler_type)",
        data=subset,
    ).fit(cov_type="HC3")

    coef = model.params["num_escalations"]
    pval = model.pvalues["num_escalations"]
    ci = model.conf_int().loc["num_escalations"]

    print(f"\n{period}")
    print(f"  Hours added per escalation event:  {coef:.2f}h")
    print(f"  95% CI:                           [{ci[0]:.2f}h, {ci[1]:.2f}h]")
    print(f"  p-value:                          {pval:.4f}")

Baseline
  Hours added per escalation event:  -1.29h
  95% CI:                           [-12.71h, 10.13h]
  p-value:                          0.8245

Both Changes
  Hours added per escalation event:  0.24h
  95% CI:                           [-3.14h, 3.63h]
  p-value:                          0.8891
Code 
cases[cases["week_number"] > 13]["sla_target_hours"].sum()
292344
Code 
res_hours = round(cases[cases["week_number"] > 13]["resolution_hours"].sum(), 2)
res_hours
258842.21
Code 
sla_new = (0.7 / res_hours) * (res_hours + (1700 * 3.95))
round(sla_new, 2)
0.72