Finding leading indicators of mobile churn with PCMCI and causal effect estimation

Churn management is a core problem for any telecom or subscription business, and it matters more as budget mobile carriers (MVNOs) gain subscribers. The standard tool is a customer-level churn score: it says who is likely to leave and feeds retention campaigns. It cannot say why company-wide churn rose this quarter, because that number moves with market conditions, competitor actions, service changes and customer experience — none of them a property of one customer. Ranking indicators by their correlation with churn is the obvious next step, and it misleads in predictable ways: an autocorrelated series inherits the correlation of a real cause, two indicators driven by the same market move together, retention spend rises because churn rose, and number-porting volume is churn that has already happened.

In this post we walk through how leading factors were separated from coincident ones. About 400 candidate variables were structured as aligned time series; PCMCI, a time-series causal discovery algorithm run with tigramite and partial-correlation tests, found the lagged parents of churn; the structure was cross-checked with the PC algorithm and score-based Bayesian-network search (Max-Min Hill-Climbing and Tabu search); and propensity score matching, IPTW and the G-formula estimated how much each driver moves churn. The result was a short list of levers with lead times and effect sizes that business teams used to set policy and marketing strategy.

Solution overview

The work is an offline analysis in Python rather than a service: a weekly panel of candidate indicators, a discovery stage that proposes lagged causes of churn, an estimation stage that sizes them, and a hand-over that turns them into priorities. The research behind it had also produced a small web application on Google Cloud Run, where a user can upload a data set, or use S&P 500 prices as an example, and run the same kind of time-series discovery and effect estimation.

The numbered steps in the diagram:

1. External market data, competitor indicators, internal business metrics and customer-experience data — about 400 variables — are collected with pandas.
2. Series of different frequencies are aligned to one weekly grid. PCMCI assumes relationships that are stable over the period, so trends and drift are dealt with before any test runs.
3. PCMCI in tigramite, with partial-correlation (ParCorr) tests, first selects a small set of lagged conditions for each variable and then tests every lagged link conditioning on the parents of both ends.
4. The PC algorithm and score-based Bayesian-network search — Max-Min Hill-Climbing and Tabu search — are run on the same variables; drivers the methods agree on carry more weight.
5. Propensity score matching, inverse probability of treatment weighting and the G-formula estimate how much each driver moves churn, adjusting for the other causes found in discovery.
6. Drivers are ranked by effect size and lead time and handed to business teams, who set policy and marketing strategies against them.

Technology stack

Step 1: Structure the candidates as aligned time series

The question was about company-level churn, so the unit of analysis is time, not the customer. About 400 potential churn-related variables were collected from four kinds of source: external market data, competitor indicators, internal business metrics and customer-experience data. They arrive at different frequencies and with different gaps, so the first job is a single panel — one row per week, one column per indicator — with a consistent rule for turning daily series into weekly ones and for carrying slower series forward.

Causal discovery on time series also assumes that relationships are stable over the period analysed. A series with a strong trend correlates with every other trending series, so long-run drift has to be removed before any test is run. In the sketch below, and in the live model, that is a linear detrend followed by standardization.

import numpy as np
import pandas as pd

def weekly_panel(sources, start, end):
    """sources: {name: pd.Series with a DatetimeIndex}, at mixed frequencies."""
    cols = {}
    for name, s in sources.items():
        s = s.loc[start:end]
        if pd.infer_freq(s.index) in ("D", "B"):
            cols[name] = s.resample("W-MON").mean()    # daily -> weekly mean
        else:
            cols[name] = s.resample("W-MON").ffill()   # monthly -> carried forward
    panel = pd.DataFrame(cols).interpolate(limit=2)
    return panel.dropna(axis=1, thresh=int(0.95 * len(panel))).dropna()

def make_stationary(panel):
    """Remove each series' linear trend, then standardize."""
    t = np.arange(len(panel))
    detrended = panel.apply(lambda s: s - np.polyval(np.polyfit(t, s, 1), t))
    return (detrended - detrended.mean()) / detrended.std()

Simplified. Real sources need their own alignment rules; a linear detrend is the simplest way to remove shared drift.

Step 2: Discover lagged parents with PCMCI

PCMCI asks, for every pair of variables and every lag up to τ_max, whether X at t − τ still carries information about Y at t once the right things are conditioned on. It works in two phases, which is what makes it workable with hundreds of variables and a few years of weekly history.

• PC1 — condition selection. For each variable, start with all lagged candidates and drop those that become independent of it when the strongest other candidates are conditioned on, growing the conditioning set step by step. A liberal threshold (pc_alpha) keeps this phase from discarding true parents; its job is a small, relevant conditioning set, not the final answer.
• MCI — momentary conditional independence. Re-test each candidate link conditioning on the parents of the effect and on the parents of the cause, shifted by the lag. Conditioning on the cause's own past removes links that only reflect autocorrelation; conditioning on the effect's other parents removes links that only reflect a common driver.

X^i_(t−τ) → X^j_t   kept iff   X^i_(t−τ) ⊥̸ X^j_t  |  P̂(X^j_t) \ {X^i_(t−τ)}  ∪  P̂(X^i_(t−τ))

ParCorr:   r = corr( X^i_(t−τ) − β̂_X·Z ,  X^j_t − β̂_Y·Z ),   Z = the conditioning set,  β̂ by least squares

The conditional-independence test is ParCorr, tigramite's partial-correlation test: regress both variables on the conditioning set and correlate the residuals. It assumes linear dependence and roughly Gaussian noise — a reasonable first approximation for aggregated weekly indicators — and is fast enough to run the thousands of tests that hundreds of variables and several lags require.

from tigramite import data_processing as pp
from tigramite.independence_tests.parcorr import ParCorr
from tigramite.pcmci import PCMCI

def churn_parents(panel, target="churn", tau_max=4, pc_alpha=0.2, alpha=0.01):
    names = list(panel.columns)
    df = pp.DataFrame(panel.to_numpy(), var_names=names)
    pcmci = PCMCI(dataframe=df, cond_ind_test=ParCorr(significance="analytic"))

    # Phase 1 (PC1): a small conditioning set per variable, liberal threshold
    parents = pcmci.run_pc_stable(tau_min=1, tau_max=tau_max, pc_alpha=pc_alpha)
    # Phase 2 (MCI): test every lagged link given the parents of both ends
    res = pcmci.run_mci(tau_min=1, tau_max=tau_max, parents=parents)
    q = pcmci.get_corrected_pvalues(p_matrix=res["p_matrix"], tau_max=tau_max,
                                    fdr_method="fdr_bh")

    j = names.index(target)
    leads = [(names[i], tau, res["val_matrix"][i, j, tau], q[i, j, tau])
             for i in range(len(names)) for tau in range(1, tau_max + 1)
             if q[i, j, tau] <= alpha]
    return sorted(leads, key=lambda r: -abs(r[2]))  # (driver, lag, partial r, q)

Simplified. Illustrated with tigramite; thresholds are illustrative, and with hundreds of variables the limits on conditioning-set size matter for run time.

Step 3: Cross-check with constraint- and score-based search

One algorithm's graph is a hypothesis. The same variables were also analysed with the PC algorithm — the cross-sectional ancestor of PCMCI, which starts from a fully connected undirected graph, removes edges by conditional-independence tests and then orients what remains — and with score-based search for a Bayesian network, which treats structure learning as optimization: search the space of directed acyclic graphs for the one with the best score.

Checking every possible graph is impossible at this size, so two heuristics were used. Max-Min Hill-Climbing first restricts each variable's candidate parents and children with a constraint-based step, then hill-climbs within that restricted space to maximize the BDeu score. Tabu search is hill-climbing that remembers recently removed or reversed edges and will not undo those moves for a while, so the search can walk out of local optima instead of cycling.

Why compare algorithms? Constraint-based and score-based methods fail in different ways — one through test thresholds, the other through search and scoring. A driver that survives both is much less likely to be an artefact of either.

Step 4: Estimate how much each driver moves churn

A graph says which arrows exist, not how large they are, and business teams needed sizes to compare levers. With candidate causes identified, effects were estimated with three methods from the potential-outcomes toolkit. Each mimics a randomized experiment in observational data, given a sufficient adjustment set — here, the other causes found in discovery, which is why discovery comes first.

• Propensity score matching pairs each treated period with an untreated one that had a similar probability of treatment given the confounders, and compares outcomes within pairs.
• Inverse probability of treatment weighting reweights periods by the inverse probability of the treatment they received, creating a pseudo-population in which treatment no longer depends on the confounders.
• The G-formula (conditional outcome modelling) fits a model of churn given treatment and confounders, then averages its predictions with treatment set on and off for every period.

G-formula   E[Y^a] = Σ_l E[Y | A = a, L = l] · P(L = l)
IPTW        ATE = E[ A·Y / e(L) ] − E[ (1 − A)·Y / (1 − e(L)) ],      e(L) = P(A = 1 | L)
PSM         match treated and untreated periods on e(L); compare outcomes

import numpy as np
from sklearn.linear_model import LinearRegression, LogisticRegression

def effects(df, treat, outcome, confounders):
    """treat: 0/1 per week; confounders: the other discovered parents, lagged."""
    A = df[treat].to_numpy()
    Y = df[outcome].to_numpy()
    L = df[confounders].to_numpy()

    # propensity score e(L) = P(A = 1 | L), clipped to avoid extreme weights
    e = LogisticRegression(max_iter=1000).fit(L, A).predict_proba(L)[:, 1]
    e = np.clip(e, 0.05, 0.95)

    # IPTW: reweight weeks by the inverse probability of what they received
    iptw = np.mean(A * Y / e) - np.mean((1 - A) * Y / (1 - e))

    # G-formula: outcome model, then average predictions with A set to 1 and 0
    om = LinearRegression().fit(np.column_stack([A, L]), Y)
    on, off = np.ones_like(A), np.zeros_like(A)
    g = np.mean(om.predict(np.column_stack([on, L])) -
                om.predict(np.column_stack([off, L])))

    # PSM: nearest untreated week on the propensity score, per treated week
    t, c = np.flatnonzero(A == 1), np.flatnonzero(A == 0)
    match = c[np.abs(e[c][None, :] - e[t][:, None]).argmin(axis=1)]
    psm = np.mean(Y[t] - Y[match])                 # effect on the treated

    return {"iptw": iptw, "g_formula": g, "psm_att": psm}

Simplified. For illustration a driver is dichotomized into high and low weeks; L holds the other discovered parents at their lags.

Matching, weighting and outcome modelling rest on different modelling assumptions — a propensity model in the first two, an outcome model in the third — so comparing them is also a robustness check: an effect with the same sign and a similar size under all three is one a business team can act on with more confidence.

Step 5: Turn drivers into levers

The output for business teams was not a graph but a short list: which indicator, how many weeks ahead of churn it moves, in which direction and by how much. Lead time decides what can be done with a driver — a factor that moves churn three weeks later leaves room to respond before the churn happens. Whether the company controls the driver decides how: internal levers can be changed directly, while external ones such as competitor activity can only be watched and prepared for.

On this basis, business teams set policy and marketing strategies for proactive churn management, concentrating limited retention resources where the analysis said they would matter. Because each driver leads churn by a known lag, the same parents also give a basis for forecasting company-level churn a few weeks ahead — one of the original goals of the project.

Feedback is the normal case. Retention spend goes up when churn goes up, so in a correlation table it looks like a cause of churn. Treated as a lagged, conditioned relationship, spend that seemed to raise churn turns out to lower it once its timing is accounted for. The live model below plants exactly this loop.

Step 6: Make discovery reusable as a web application

The methods used here came out of research on time-series causal discovery that had already been turned into a minimum viable web application on Google Cloud Run, so that the analysis was not tied to one data set or one analyst. A user uploads a data set — or uses S&P 500 stock prices as an example, to ask which companies' prices lead a chosen company's — selects the time index, the target and the candidate variables, sets the minimum and maximum lag and the test strength, and gets back the causal graph together with each candidate's causal effect and direct effect on the target.

The two effects answer different questions. The direct effect is the strength of the link itself; the total causal effect also counts paths through other variables. A driver that moves churn mostly by first moving a second indicator has a small direct effect and a large total one — and acting on it still works.

from tigramite import data_processing as pp
from tigramite.independence_tests.parcorr import ParCorr
from tigramite.models import LinearMediation
from tigramite.pcmci import PCMCI

def analyse(data, names, target, candidates, tau_min, tau_max, alpha):
    """What the MVP returns once the user has chosen target, candidates and lags."""
    df = pp.DataFrame(data, var_names=names)
    pcmci = PCMCI(dataframe=df, cond_ind_test=ParCorr())
    res = pcmci.run_pcmci(tau_min=tau_min, tau_max=tau_max, alpha_level=alpha)
    parents = pcmci.return_parents_dict(graph=res["graph"],
                                        val_matrix=res["val_matrix"])

    med = LinearMediation(dataframe=df)
    med.fit_model(all_parents=parents, tau_max=tau_max)
    j = names.index(target)
    effects = [{"cause": c, "lag": tau,
                "total": med.get_ce(i=names.index(c), tau=-tau, j=j),
                "direct": med.get_coeff(i=names.index(c), tau=-tau, j=j)}
               for c in candidates
               for tau in range(max(tau_min, 1), tau_max + 1)]
    return res["graph"], effects

Simplified. Illustrated with tigramite's linear mediation model; upload handling, validation and plotting are omitted.

Results

The analysis supported policy and marketing strategies for proactive churn management and contributed to reducing the churn-rate gap against the leading competitor. More durably, it expanded churn management from customer-level prediction into company-level leading-factor management.

The approach also travelled beyond churn: the same time-series causal-discovery application later supported the search for process conditions that move manufacturing yield.

Lessons learned

• Ask about time, not just association. A cause must precede its effect. Testing lagged links while conditioning on both variables' pasts removed the two classic false positives of correlation analysis — autocorrelation and common drivers.
• Select liberally, test strictly. PCMCI's split — a permissive phase to build conditioning sets, a strict phase to decide links — keeps power with hundreds of candidates and limited history.
• Compare algorithms, trust agreement. Constraint-based and score-based searches fail differently; drivers they agreed on carried more weight.
• Hand over sizes and lead times, not arrows. Effect estimates turned a graph into a list of levers that teams could compare and budget against.

Conclusion

Company-level churn is driven by markets, competitors, service and experience, and a customer-level score cannot see any of it. Structuring hundreds of indicators as time series, discovering lagged causes with PCMCI, cross-checking with other structure-learning methods and sizing effects with matching, weighting and outcome models turned a correlation table into a short list of leading factors that business teams acted on.

The pattern — candidate series, lagged discovery, adjustment sets taken from the graph, effect estimates for decisions — applies wherever an aggregate rate moves with many indicators and the question is which of them to act on.

Limitations

• Observational causal discovery rests on assumptions — no unmeasured common causes, stable relationships, the right lag window — that can be argued but not proven; effect estimates are evidence for a decision, not a guarantee of its outcome.
• Short histories limit what can be found: with a few years of weekly data, weak or long-lag drivers are missed and a few false links survive even strict thresholds.
• Partial-correlation tests capture linear dependence; non-linear or threshold effects need other conditional-independence tests.
• The live model uses nine generated indicators, linear detrending, partial-correlation tests and a linear outcome model on the discovered parents; the project's ~400 variables, its other discovery algorithms and its matching and weighting estimators are not reproduced.

About the demo and confidentiality

The indicators, their names, the planted relationships and all numbers in the embedded model are invented. No market, competitor, operational or customer data from any operator appears here, the real variable set and the drivers found are not described, and code is simplified and written for illustration.