Correlation, Causation and the Missing Variable

import numpy as np
import pandas as pd

rng = np.random.default_rng(3)
n = 1000

# A pure confounder: company size drives BOTH support tickets and
# revenue. Neither causes the other.
size = rng.lognormal(3.0, 0.8, n)
tickets = size * 0.9 + rng.normal(0, 2, n)
revenue = size * 140 + rng.normal(0, 200, n)

frame = pd.DataFrame({'size': size, 'tickets': tickets, 'revenue': revenue})
print(frame[['tickets', 'revenue']].corr().round(3))
print()
print(f'tickets vs revenue: r = {frame["tickets"].corr(frame["revenue"]):.3f}')
print('"Customers who raise more tickets spend more. Should we')
print('encourage ticket-raising?"')

import numpy as np
import pandas as pd

rng = np.random.default_rng(3)
n = 1000
size = rng.lognormal(3.0, 0.8, n)
tickets = size * 0.9 + rng.normal(0, 2, n)
revenue = size * 140 + rng.normal(0, 200, n)
frame = pd.DataFrame({'size': size, 'tickets': tickets, 'revenue': revenue})

# Within a narrow band of company size, is there still a relationship?
frame['band'] = pd.qcut(frame['size'], 5, labels=[f'q{i}' for i in range(1, 6)])
within = frame.groupby('band', observed=True).apply(
    lambda g: g['tickets'].corr(g['revenue']), include_groups=False)

print(f'overall correlation        {frame["tickets"].corr(frame["revenue"]):+.3f}')
print('within each size quintile:')
print(within.round(3).to_string())
print()
print(f'mean within-band correlation {within.mean():+.3f}')
print('The relationship was entirely company size.')

import pandas as pd

claims = pd.DataFrame([
    {'claim': 'users who use feature X retain better',
     'reverse': 'engaged users find X; X does not create engagement',
     'confounder': 'tenure - older accounts use more features and retain more'},
    {'claim': 'customers who call support churn more',
     'reverse': 'people about to churn call support',
     'confounder': 'product problems cause both'},
    {'claim': 'the email campaign lifted sales',
     'reverse': 'intent to buy drives opening the email',
     'confounder': 'the campaign ran in the week of a holiday'},
    {'claim': 'teams using our tool ship faster',
     'reverse': 'fast teams adopt tools',
     'confounder': 'engineering maturity causes both'},
    {'claim': 'the redesign increased sign-ups',
     'reverse': 'none plausible',
     'confounder': 'a press mention in the same week'},
])
for _, row in claims.iterrows():
    print(f'claim      {row["claim"]}')
    print(f'  reverse  {row["reverse"]}')
    print(f'  confound {row["confounder"]}')
    print()

import numpy as np
import pandas as pd

rng = np.random.default_rng(11)
n = 5000

# Two genuinely INDEPENDENT qualities of a candidate.
skill = rng.normal(0, 1, n)
interview = rng.normal(0, 1, n)
print(f'in the whole population, r = {np.corrcoef(skill, interview)[0, 1]:+.3f}')

# Hire anyone who is good at EITHER. Among the hired, the two become
# negatively correlated - with nothing having caused anything.
hired = (skill + interview) > 1.2
print(f'among those hired,      r = '
      f'{np.corrcoef(skill[hired], interview[hired])[0, 1]:+.3f}')
print(f'({hired.sum()} of {n} hired)')
print()
print('"Our best interviewers are the weakest engineers" is an artefact')
print('of the hiring rule, not a fact about people.')

import numpy as np
import pandas as pd

rng = np.random.default_rng(17)

# Two regions. The feature ships in North only, in week 13.
weeks = np.arange(1, 27)
rows = []
for region, base, effect in [('North', 100, 12), ('South', 85, 0)]:
    for week in weeks:
        seasonal = 8 * np.sin(week / 4)          # affects BOTH regions
        treated = effect if (region == 'North' and week >= 13) else 0
        rows.append({'region': region, 'week': week,
                     'value': base + seasonal + treated + rng.normal(0, 3)})
data = pd.DataFrame(rows)
data['period'] = np.where(data['week'] >= 13, 'after', 'before')

cell = data.pivot_table(index='region', columns='period', values='value', aggfunc='mean')
cell['change'] = cell['after'] - cell['before']
print(cell.round(1).to_string())
print()

naive = cell.loc['North', 'change']
did = cell.loc['North', 'change'] - cell.loc['South', 'change']
print(f'before-and-after in North only: {naive:+.1f}')
print(f'difference-in-differences:      {did:+.1f}')
print(f'the true effect is:             +12.0')
print()
print('The naive number includes the seasonal swing that hit both')
print('regions. Subtracting the control removes it.')

import numpy as np
import pandas as pd

rng = np.random.default_rng(23)
n = 3000

users = pd.DataFrame({
    'user_id': np.arange(n),
    'signup_week': rng.integers(1, 20, n),
})
users['active_weeks'] = rng.poisson(6, n) + 1
users['retained'] = users['active_weeks'] >= 8

# The broken measure: feature use counted over the user's WHOLE life.
users['feature_uses_lifetime'] = rng.poisson(users['active_weeks'] * 0.8)

# The honest measure: feature use in the FIRST TWO WEEKS only.
users['feature_uses_first_2w'] = rng.poisson(1.6, n)

for column in ['feature_uses_lifetime', 'feature_uses_first_2w']:
    used = users[column] > 0
    rate_used = users.loc[used, 'retained'].mean()
    rate_not = users.loc[~used, 'retained'].mean()
    print(f'{column}')
    print(f'  used     {rate_used:.1%} retained  (n={used.sum():,})')
    print(f'  not used {rate_not:.1%} retained  (n={(~used).sum():,})')
    print(f'  gap      {rate_used - rate_not:+.1%}')
    print()
print('The lifetime measure is circular: retained users had more weeks')
print('in which to use the feature. The two-week measure is not.')

import pandas as pd

print(pd.DataFrame([
    {'instead of': 'X drives retention',
     'say': 'X use is associated with higher retention; we have not tested whether it causes it'},
    {'instead of': 'the campaign increased sales 12%',
     'say': 'sales were 12% higher in campaign weeks; a holiday fell in the same window'},
    {'instead of': 'support calls cause churn',
     'say': 'churned accounts called support more in their final month'},
    {'instead of': 'we should make everyone use X',
     'say': 'an experiment on X would cost two weeks and would answer this'},
]).to_string(index=False))

import numpy as np
import pandas as pd

rng = np.random.default_rng(seed=167)

# A product team believes: "users who join a team workspace retain
# better, so we should push everyone into workspaces."
#
# The data is generated so we KNOW the truth: workspace membership has
# a small real effect (+4 points), and company size drives both
# workspace use and retention much more strongly.

n = 6000
company_size = rng.lognormal(2.4, 0.9, n)                  # the confounder
size_band = pd.qcut(company_size, 4, labels=['tiny', 'small', 'medium', 'large'])

# Bigger companies are far likelier to use a workspace...
workspace_probability = 0.08 + 0.72 * (company_size / (company_size + 12))
in_workspace = rng.random(n) < workspace_probability

# ...and bigger companies retain much better, independently.
retention_probability = (0.22 + 0.55 * (company_size / (company_size + 10))
                         + 0.04 * in_workspace)            # the TRUE effect
retained = rng.random(n) < retention_probability

users = pd.DataFrame({
    'company_size': company_size, 'band': size_band,
    'workspace': in_workspace, 'retained': retained,
})

# ---- 1. The naive comparison, which is what gets presented. --------
naive = users.groupby('workspace')['retained'].agg(['mean', 'size'])
gap = naive.loc[True, 'mean'] - naive.loc[False, 'mean']
print('naive comparison:')
print(naive.assign(mean=lambda d: (d['mean'] * 100).round(1)).to_string())
print(f'gap {gap:+.1%} - "workspaces lift retention by {gap * 100:.0f} points"')
print(f'the TRUE effect built into this data is +4.0 points')
print()

# ---- 2. Name the confounder, then stratify. ------------------------
within = users.groupby(['band', 'workspace'], observed=True)['retained'].agg(
    ['mean', 'size']).unstack()
within.columns = ['no workspace', 'workspace', 'n no', 'n yes']
within['gap_points'] = ((within['workspace'] - within['no workspace']) * 100).round(1)
print('within each company-size band:')
print(within.assign(**{c: (within[c] * 100).round(1)
                       for c in ['no workspace', 'workspace']}).to_string())
print()
weighted = np.average(within['gap_points'], weights=within['n no'] + within['n yes'])
print(f'size-weighted gap within bands: {weighted:+.1f} points')
print(f'naive gap:                      {gap * 100:+.1f} points')
print(f'true effect:                    +4.0 points')
print()
print('Most of the headline gap was company size. The within-band')
print('estimate is close to the truth.')
print()

# ---- 3. Reverse causation: is it plausible here? -------------------
print('reverse causation check:')
print('  Could retaining cause workspace use rather than the reverse?')
print('  Yes - a user who stays has more opportunity to be invited.')
print('  The stratified estimate does not rule that out.')
print()

# ---- 4. Selection: who is even in this data? -----------------------
print('selection check:')
print(f'  {len(users):,} users, {users["workspace"].mean():.1%} in a workspace')
print('  Users who left before day 7 are not in the retention table at')
print('  all, and they are disproportionately small companies.')
print()

# ---- 5. What would settle it, with a cost. -------------------------
BASELINE = float(users.loc[~users['workspace'], 'retained'].mean())
def sample_size(baseline, absolute_lift):
    z_alpha, z_beta = 1.959964, 0.841621
    p1, p2 = baseline, baseline + absolute_lift
    p_bar = (p1 + p2) / 2
    return int(np.ceil(((z_alpha * np.sqrt(2 * p_bar * (1 - p_bar))
                         + z_beta * np.sqrt(p1 * (1 - p1) + p2 * (1 - p2))) ** 2)
                       / (p2 - p1) ** 2))

per_group = sample_size(BASELINE, 0.04)
print('what would settle it:')
print(f'  Randomly prompt half of new users to create a workspace.')
print(f'  To detect a 4-point lift from a {BASELINE:.0%} baseline needs')
print(f'  {per_group:,} users per group, {2 * per_group:,} in total.')
print()

# ---- 6. The sentence, written honestly. ----------------------------
print('Workspace members retain '
      f'{gap * 100:.0f} points better than non-members '
      f'({naive.loc[True, "mean"]:.0%} against {naive.loc[False, "mean"]:.0%}), '
      'but company size explains most of that: within size bands the gap '
      f'falls to {weighted:.0f} points. We have not ruled out that staying '
      'causes workspace membership rather than the reverse. A randomised '
      f'prompt to {2 * per_group:,} new users would answer it in a few weeks.')