生存時間分析をSQLで近似する方法｜顧客が離脱するまでの期間を測る

S(t) = P(T > t) ＝ 時点 t より後まで「生存」している確率

-- orders テーブル
-- order_id      : 注文ID
-- customer_id   : 顧客ID
-- order_date    : 注文日
-- status        : 'completed' / 'cancelled'

-- STEP1: 生存時間とイベントフラグの計算

WITH customer_orders AS (
    SELECT
        customer_id,
        MIN(order_date)  AS first_order_date,
        MAX(order_date)  AS last_order_date,
        COUNT(order_id)  AS purchase_count
    FROM orders
    WHERE status = 'completed'
    GROUP BY customer_id
),
survival_data AS (
    SELECT
        customer_id,
        first_order_date,
        last_order_date,
        purchase_count,

        -- 基準日
        DATE '2024-12-31'  AS reference_date,

        -- 最終購買からの経過日数
        DATE_DIFF('day', last_order_date, DATE '2024-12-31')
            AS days_since_last_order,

        -- イベント（離脱）の発生判定
        -- 180日以上沈黙 = 離脱イベント発生
        CASE
            WHEN DATE_DIFF('day', last_order_date, DATE '2024-12-31') >= 180
            THEN 1 ELSE 0
        END  AS event_occurred,

        -- 生存時間（duration）
        -- イベント発生の場合：初回購買から離脱判定日まで
        -- 打ち切りの場合：初回購買から基準日まで
        CASE
            WHEN DATE_DIFF('day', last_order_date, DATE '2024-12-31') >= 180
            -- 離脱した場合、「最終購買日 + 180日」が離脱判定日
            THEN DATE_DIFF('day', first_order_date,
                     DATE_ADD('day', 180, last_order_date))
            -- 打ち切りの場合、基準日までの全期間
            ELSE DATE_DIFF('day', first_order_date, DATE '2024-12-31')
        END  AS duration_days
    FROM customer_orders
)
SELECT
    customer_id,
    first_order_date,
    last_order_date,
    purchase_count,
    days_since_last_order,
    event_occurred,
    duration_days,
    CASE WHEN event_occurred = 1 THEN '離脱' ELSE '打ち切り（継続中）' END  AS status_label
FROM survival_data
ORDER BY duration_days ASC;

各時点 t ごとに：
  n_t = その時点で「生存中」だった顧客数（リスク集合）
  d_t = その時点で「離脱」した顧客数
  S(t) = S(t-1) × (1 - d_t / n_t)

-- STEP2: Kaplan-Meier 推定量の計算

WITH customer_orders AS (
    SELECT
        customer_id,
        MIN(order_date) AS first_order_date,
        MAX(order_date) AS last_order_date
    FROM orders
    WHERE status = 'completed'
    GROUP BY customer_id
),
survival_data AS (
    SELECT
        customer_id,
        CASE
            WHEN DATE_DIFF('day', last_order_date, DATE '2024-12-31') >= 180 THEN 1
            ELSE 0
        END  AS event_occurred,
        CASE
            WHEN DATE_DIFF('day', last_order_date, DATE '2024-12-31') >= 180
            THEN DATE_DIFF('day', first_order_date, DATE_ADD('day', 180, last_order_date))
            ELSE DATE_DIFF('day', first_order_date, DATE '2024-12-31')
        END  AS duration_days
    FROM customer_orders
),
-- 時点ごとのイベント数・打ち切り数を集計
-- 30日単位のビンに丸める（日次だと行数が多すぎるため）
event_table AS (
    SELECT
        -- 30日単位のビン
        FLOOR(duration_days / 30) * 30  AS time_bin,
        SUM(event_occurred)              AS events,       -- その期間の離脱数
        COUNT(*) - SUM(event_occurred)   AS censored      -- その期間の打ち切り数
    FROM survival_data
    GROUP BY FLOOR(duration_days / 30) * 30
),
-- リスク集合の計算（各時点での「生存中」の人数）
risk_set AS (
    SELECT
        time_bin,
        events,
        censored,
        -- リスク集合 = 全体 - その時点より前に離脱・打ち切りになった人数
        SUM(events + censored) OVER (
            ORDER BY time_bin DESC
            ROWS BETWEEN UNBOUNDED PRECEDING AND CURRENT ROW
        )  AS at_risk
    FROM event_table
),
-- Kaplan-Meier 推定量の計算
km_estimate AS (
    SELECT
        time_bin,
        at_risk,
        events,
        censored,
        -- 各時点の生存確率の積（KM推定の核心）
        EXP(SUM(LN(CASE WHEN at_risk > 0
                   THEN 1.0 - events * 1.0 / at_risk
                   ELSE 1.0 END))
            OVER (ORDER BY time_bin ASC
                  ROWS BETWEEN UNBOUNDED PRECEDING AND CURRENT ROW)
        )  AS survival_prob,
        -- Greenwood の公式による95%信頼区間の近似
        -- Var(log S(t)) ≈ Σ d_i / [n_i(n_i - d_i)]
        SUM(
            CASE WHEN at_risk > 1 AND events > 0
                 THEN events * 1.0 / (at_risk * (at_risk - events))
                 ELSE 0 END
        ) OVER (ORDER BY time_bin ASC
                ROWS BETWEEN UNBOUNDED PRECEDING AND CURRENT ROW)  AS greenwood_var
    FROM risk_set
)
SELECT
    time_bin                                                AS days_since_first_order,
    at_risk,
    events                                                  AS churn_events,
    censored,
    ROUND(survival_prob, 4)                                 AS survival_rate,
    ROUND((1 - survival_prob) * 100, 1)                    AS churn_rate_pct,

    -- 95%信頼区間（Greenwood公式 + 対数変換）
    ROUND(GREATEST(0, survival_prob * EXP(-1.96 * SQRT(greenwood_var))), 4)  AS ci_lower,
    ROUND(LEAST(1,   survival_prob * EXP( 1.96 * SQRT(greenwood_var))), 4)   AS ci_upper
FROM km_estimate
ORDER BY time_bin;

-- STEP3: チャネル別 Kaplan-Meier 比較

WITH customer_base AS (
    SELECT
        o.customer_id,
        c.channel,
        MIN(o.order_date)  AS first_order_date,
        MAX(o.order_date)  AS last_order_date
    FROM orders     o
    JOIN customers  c  ON o.customer_id = c.customer_id
    WHERE o.status = 'completed'
    GROUP BY o.customer_id, c.channel
),
survival_data AS (
    SELECT
        customer_id,
        channel,
        CASE
            WHEN DATE_DIFF('day', last_order_date, DATE '2024-12-31') >= 180 THEN 1
            ELSE 0
        END  AS event_occurred,
        CASE
            WHEN DATE_DIFF('day', last_order_date, DATE '2024-12-31') >= 180
            THEN DATE_DIFF('day', first_order_date, DATE_ADD('day', 180, last_order_date))
            ELSE DATE_DIFF('day', first_order_date, DATE '2024-12-31')
        END  AS duration_days
    FROM customer_base
),
event_table AS (
    SELECT
        channel,
        FLOOR(duration_days / 90) * 90  AS time_bin,  -- 90日単位
        SUM(event_occurred)              AS events,
        COUNT(*) - SUM(event_occurred)   AS censored
    FROM survival_data
    GROUP BY channel, FLOOR(duration_days / 90) * 90
),
risk_set AS (
    SELECT
        channel,
        time_bin,
        events,
        SUM(events + censored) OVER (
            PARTITION BY channel
            ORDER BY time_bin DESC
            ROWS BETWEEN UNBOUNDED PRECEDING AND CURRENT ROW
        )  AS at_risk
    FROM event_table
)
SELECT
    channel,
    time_bin  AS days,
    at_risk,
    events,
    ROUND(
        EXP(SUM(LN(CASE WHEN at_risk > 0
                   THEN 1.0 - events * 1.0 / at_risk
                   ELSE 1.0 END))
            OVER (PARTITION BY channel
                  ORDER BY time_bin ASC
                  ROWS BETWEEN UNBOUNDED PRECEDING AND CURRENT ROW))
    , 4)  AS survival_rate
FROM risk_set
ORDER BY channel, time_bin;

-- STEP4: ハザード率の計算（どの時点で最も離脱リスクが高いか）

WITH km_data AS (
    -- STEP2の集計済みテーブルを使う（ここでは簡略化）
    SELECT
        time_bin,
        at_risk,
        events,
        ROUND(
            EXP(SUM(LN(CASE WHEN at_risk > 0
                       THEN 1.0 - events * 1.0 / at_risk
                       ELSE 1.0 END))
                OVER (ORDER BY time_bin ASC
                      ROWS BETWEEN UNBOUNDED PRECEDING AND CURRENT ROW))
        , 4)  AS survival_rate
    FROM (
        -- STEP2のrisk_set相当の集計
        SELECT
            FLOOR(duration_days / 30) * 30  AS time_bin,
            SUM(events + censored) OVER (ORDER BY FLOOR(duration_days / 30) * 30 DESC
                                         ROWS BETWEEN UNBOUNDED PRECEDING AND CURRENT ROW) AS at_risk,
            SUM(event_occurred)  AS events,
            COUNT(*) - SUM(event_occurred)  AS censored
        FROM (
            SELECT
                CASE WHEN DATE_DIFF('day', MAX(order_date), DATE '2024-12-31') >= 180 THEN 1 ELSE 0 END AS event_occurred,
                CASE WHEN DATE_DIFF('day', MAX(order_date), DATE '2024-12-31') >= 180
                     THEN DATE_DIFF('day', MIN(order_date), DATE_ADD('day', 180, MAX(order_date)))
                     ELSE DATE_DIFF('day', MIN(order_date), DATE '2024-12-31') END AS duration_days
            FROM orders WHERE status='completed' GROUP BY customer_id
        )
        GROUP BY FLOOR(duration_days / 30) * 30
    )
)
SELECT
    time_bin  AS days,
    at_risk,
    events,
    survival_rate,
    -- ハザード率 = 各時点での離脱数 / リスク集合
    ROUND(events * 1.0 / NULLIF(at_risk, 0), 4)  AS hazard_rate,
    -- 前の期間との比較
    LAG(events * 1.0 / NULLIF(at_risk, 0)) OVER (ORDER BY time_bin)  AS prev_hazard
FROM km_data
ORDER BY time_bin;