{"spec_id":"survival-kaplan-meier","library":"makie","language":"julia","code":"# anyplot.ai\n# survival-kaplan-meier: Kaplan-Meier Survival Plot\n# Library: makie 0.21.9 | Julia 1.11.9\n# Quality: 92/100 | Created: 2026-09-09\n\nusing CairoMakie\nusing Colors\nusing Random\n\nRandom.seed!(42)\n\n# Theme tokens\nconst THEME       = get(ENV, \"ANYPLOT_THEME\", \"light\")\nconst PAGE_BG     = THEME == \"light\" ? colorant\"#FAF8F1\" : colorant\"#1A1A17\"\nconst ELEVATED_BG = THEME == \"light\" ? colorant\"#FFFDF6\" : colorant\"#242420\"\nconst INK         = THEME == \"light\" ? colorant\"#1A1A17\" : colorant\"#F0EFE8\"\nconst INK_SOFT    = THEME == \"light\" ? colorant\"#4A4A44\" : colorant\"#B8B7B0\"\n\nconst IMPRINT_PALETTE = [\n    colorant\"#009E73\",\n    colorant\"#C475FD\",\n    colorant\"#4467A3\",\n    colorant\"#BD8233\",\n    colorant\"#AE3030\",\n    colorant\"#2ABCCD\",\n    colorant\"#954477\",\n    colorant\"#99B314\",\n]\n\n# Data — subscription retention: Free vs Premium tier time-to-churn (months)\nn_free    = 220\nn_premium = 220\nstudy_end = 36.0\n\nchurn_free     = 14.0 .* (-log.(rand(n_free))) .^ (1 / 1.3)\nchurn_premium  = 26.0 .* (-log.(rand(n_premium))) .^ (1 / 1.3)\ncensor_free    = min.(study_end, 40.0 .* (-log.(rand(n_free))))\ncensor_premium = min.(study_end, 40.0 .* (-log.(rand(n_premium))))\n\ntime_free     = min.(churn_free, censor_free)\nevent_free    = Int.(churn_free .<= censor_free)\ntime_premium  = min.(churn_premium, censor_premium)\nevent_premium = Int.(churn_premium .<= censor_premium)\n\n# Kaplan-Meier estimator with Greenwood's formula for the standard error\nfunction kaplan_meier(time, event)\n    event_times = sort(unique(time[event .== 1]))\n    survival = Float64[]\n    stderr   = Float64[]\n\n    s = 1.0\n    greenwood_sum = 0.0\n    for et in event_times\n        n_risk = count(>=(et), time)\n        d      = count(==(et), time[event .== 1])\n        s *= (1 - d / n_risk)\n        push!(survival, s)\n        if n_risk > d\n            greenwood_sum += d / (n_risk * (n_risk - d))\n        end\n        push!(stderr, s * sqrt(greenwood_sum))\n    end\n\n    censor_times = time[event .== 0]\n    return event_times, survival, stderr, censor_times\nend\n\n# Step-function coordinates for the curve and its 95% confidence band\nfunction km_step(event_times, survival, stderr, t_start, t_end)\n    xs, ys, lo, hi = Float64[t_start], Float64[1.0], Float64[1.0], Float64[1.0]\n    for i in eachindex(event_times)\n        push!(xs, event_times[i]); push!(ys, ys[end]); push!(lo, lo[end]); push!(hi, hi[end])\n        push!(xs, event_times[i]); push!(ys, survival[i])\n        push!(lo, clamp(survival[i] - 1.96 * stderr[i], 0.0, 1.0))\n        push!(hi, clamp(survival[i] + 1.96 * stderr[i], 0.0, 1.0))\n    end\n    push!(xs, t_end); push!(ys, ys[end]); push!(lo, lo[end]); push!(hi, hi[end])\n    return xs, ys, lo, hi\nend\n\nsurvival_at(event_times, survival, t) = isempty(event_times) || t < event_times[1] ? 1.0 :\n    survival[findlast(<=(t), event_times)]\n\net_free, surv_free, se_free, cens_free = kaplan_meier(time_free, event_free)\net_prem, surv_prem, se_prem, cens_prem = kaplan_meier(time_premium, event_premium)\n\nxs_free, ys_free, lo_free, hi_free = km_step(et_free, surv_free, se_free, 0.0, study_end)\nxs_prem, ys_prem, lo_prem, hi_prem = km_step(et_prem, surv_prem, se_prem, 0.0, study_end)\n\ncens_y_free = [survival_at(et_free, surv_free, t) for t in cens_free]\ncens_y_prem = [survival_at(et_prem, surv_prem, t) for t in cens_prem]\n\nmedian_idx_free = findfirst(<=(0.5), surv_free)\nmedian_idx_prem = findfirst(<=(0.5), surv_prem)\nmedian_free = median_idx_free === nothing ? nothing : et_free[median_idx_free]\nmedian_prem = median_idx_prem === nothing ? nothing : et_prem[median_idx_prem]\n\n# Log-rank test (Mantel-Haenszel chi-square, 1 df) comparing the two groups\ncombined_event_times = sort(unique(vcat(et_free, et_prem)))\nobserved_free, expected_free, variance_sum = 0.0, 0.0, 0.0\nfor t in combined_event_times\n    n1 = count(>=(t), time_free)\n    n2 = count(>=(t), time_premium)\n    d1 = count(==(t), time_free[event_free .== 1])\n    d2 = count(==(t), time_premium[event_premium .== 1])\n    n, d = n1 + n2, d1 + d2\n    if n > 1 && d > 0\n        global observed_free += d1\n        global expected_free += d * n1 / n\n        global variance_sum += d * (n1 / n) * (n2 / n) * ((n - d) / (n - 1))\n    end\nend\nlogrank_chi2 = (observed_free - expected_free)^2 / variance_sum\n\n# Standard normal CDF (Abramowitz & Stegun 26.2.17) — avoids a SpecialFunctions dependency\nfunction normal_cdf(x)\n    z = abs(x)\n    t = 1 / (1 + 0.2316419 * z)\n    poly = t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 +\n           t * (-1.821255978 + t * 1.330274429))))\n    phi = 1 - poly * exp(-z^2 / 2) / sqrt(2π)\n    return x >= 0 ? phi : 1 - phi\nend\nlogrank_p = 2 * (1 - normal_cdf(sqrt(logrank_chi2)))\n\n# Number-at-risk table gridpoints\nrisk_times = collect(0.0:6.0:study_end)\nrisk_free  = [count(>=(t), time_free) for t in risk_times]\nrisk_prem  = [count(>=(t), time_premium) for t in risk_times]\n\n# Title (scaled fontsize — descriptive prefix pushes past the 67-char baseline)\ntitle_str  = \"Subscription Retention · survival-kaplan-meier · julia · makie · anyplot.ai\"\ntitle_size = max(14, round(Int, 20 * min(1.0, 67.0 / length(title_str))))\n\n# Plot\nfig = Figure(\n    size            = (1600, 900),\n    fontsize        = 14,\n    backgroundcolor = PAGE_BG,\n)\n\nax = Axis(\n    fig[1, 1];\n    title             = title_str,\n    titlesize         = title_size,\n    titlecolor        = INK,\n    ylabel            = \"Survival Probability\",\n    xlabelsize        = 14,\n    ylabelsize        = 14,\n    xticklabelsize    = 12,\n    yticklabelsize    = 12,\n    xlabelcolor       = INK,\n    ylabelcolor       = INK,\n    xticklabelcolor   = INK_SOFT,\n    yticklabelcolor   = INK_SOFT,\n    xtickcolor        = INK_SOFT,\n    ytickcolor        = INK_SOFT,\n    backgroundcolor   = PAGE_BG,\n    topspinevisible   = false,\n    rightspinevisible = false,\n    leftspinecolor    = INK_SOFT,\n    bottomspinecolor  = INK_SOFT,\n    xgridcolor        = RGBAf(INK.r, INK.g, INK.b, 0.12),\n    ygridcolor        = RGBAf(INK.r, INK.g, INK.b, 0.12),\n    xminorgridvisible = false,\n    yminorgridvisible = false,\n    limits            = (0.0, study_end, 0.0, 1.05),\n)\n\n# 95% confidence bands (Greenwood's formula) — drawn first, behind the curves\nband!(ax, xs_free, lo_free, hi_free;\n    color = RGBAf(IMPRINT_PALETTE[1].r, IMPRINT_PALETTE[1].g, IMPRINT_PALETTE[1].b, 0.15))\nband!(ax, xs_prem, lo_prem, hi_prem;\n    color = RGBAf(IMPRINT_PALETTE[2].r, IMPRINT_PALETTE[2].g, IMPRINT_PALETTE[2].b, 0.15))\n\n# Kaplan-Meier step curves\nlines!(ax, xs_free, ys_free; color = IMPRINT_PALETTE[1], linewidth = 3.0, label = \"Free Tier\")\nlines!(ax, xs_prem, ys_prem; color = IMPRINT_PALETTE[2], linewidth = 3.0, label = \"Premium Tier\")\n\n# Censoring marks\nscatter!(ax, cens_free, cens_y_free;\n    marker = :vline, markersize = 14, color = IMPRINT_PALETTE[1], strokewidth = 0)\nscatter!(ax, cens_prem, cens_y_prem;\n    marker = :vline, markersize = 14, color = IMPRINT_PALETTE[2], strokewidth = 0)\n\n# Median survival reference lines\nhlines!(ax, [0.5]; color = INK_SOFT, linewidth = 1.0, linestyle = :dot)\nmedian_free !== nothing && vlines!(ax, [median_free]; color = IMPRINT_PALETTE[1], linewidth = 1.0, linestyle = :dash)\nmedian_prem !== nothing && vlines!(ax, [median_prem]; color = IMPRINT_PALETTE[2], linewidth = 1.0, linestyle = :dash)\n\nmedian_free_str = median_free === nothing ? \"not reached\" : \"$(round(median_free; digits = 1)) mo\"\nmedian_prem_str = median_prem === nothing ? \"not reached\" : \"$(round(median_prem; digits = 1)) mo\"\np_str = logrank_p < 0.001 ? \"p < 0.001\" : \"p = $(round(logrank_p; digits = 3))\"\n\ntext!(ax,\n    \"Median survival — Free: $(median_free_str) · Premium: $(median_prem_str)\\nLog-rank test: $(p_str)\";\n    position = Point2f(study_end * 0.97, 0.8),\n    align    = (:right, :top),\n    color    = INK,\n    fontsize = 13,\n)\n\naxislegend(ax;\n    position        = :rt,\n    labelsize       = 12,\n    framevisible    = true,\n    framecolor      = INK_SOFT,\n    backgroundcolor = ELEVATED_BG,\n    labelcolor      = INK,\n)\n\n# Number-at-risk table\nax_risk = Axis(\n    fig[2, 1];\n    backgroundcolor     = PAGE_BG,\n    limits              = (0.0, study_end, 0.0, 2.0),\n    xlabel              = \"Time Since Signup (months)\",\n    xlabelsize          = 14,\n    xlabelcolor         = INK,\n    xticklabelsize      = 12,\n    xticklabelcolor     = INK_SOFT,\n    xtickcolor          = INK_SOFT,\n    yticks              = ([0.5, 1.5], [\"Premium\", \"Free\"]),\n    yticklabelsize      = 12,\n    yticklabelcolor     = INK_SOFT,\n    yticksvisible       = false,\n    topspinevisible     = false,\n    rightspinevisible   = false,\n    leftspinevisible    = false,\n    bottomspinecolor    = INK_SOFT,\n)\nLabel(fig[2, 1, Top()], \"Number at risk\";\n    fontsize = 12, color = INK_SOFT, halign = :left, padding = (0, 0, 4, 0))\n\nfor (i, t) in enumerate(risk_times)\n    halign = i == 1 ? :left : i == length(risk_times) ? :right : :center\n    text!(ax_risk, string(risk_free[i]);\n        position = Point2f(t, 1.5), align = (halign, :center), color = IMPRINT_PALETTE[1], fontsize = 12)\n    text!(ax_risk, string(risk_prem[i]);\n        position = Point2f(t, 0.5), align = (halign, :center), color = IMPRINT_PALETTE[2], fontsize = 12)\nend\n\nrowsize!(fig.layout, 2, Relative(0.16))\n\n# Save\nsave(\"plot-$(THEME).png\", fig; px_per_unit = 2)\n"}