Shares by country

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import plotstyle
from plotstyle import tidy

REPO = Path("..")  # the CSV sits in the repository root, above book/
VERDICTS = [
    "likely_monetized",
    "possibly_monetized",
    "likely_not_monetized",
    "unable_to_determine",
]

# keep_default_na=False: "NA" is Namibia's code, not a missing value
snapshot = pd.read_csv(
    REPO / "youtube-monetization-signals-by-country.csv", keep_default_na=False
)
snapshot.head()

checked = snapshot["channels_checked"]
assert (snapshot[VERDICTS].sum(axis=1) == checked).all()
assert checked.min() >= 30

countries = len(snapshot)
channels = int(checked.sum())
smallest, largest = int(checked.min()), int(checked.max())
rows_at_60 = int((checked == 60).sum())

# Figures quoted on the text pages of this book
assert (countries, channels) == (50, 3080)
# "one channel in five" could not be determined
assert round(snapshot["unable_to_determine"].sum() / channels, 2) == 0.20
us = snapshot.set_index("country_code").loc["US"]
assert (us["unable_to_determine"], us["channels_checked"]) == (16, 66)

print(f"{countries} countries, {channels:,} channels, "
      f"{smallest} to {largest} per country")

pooled = snapshot[VERDICTS].sum()
share_text = {
    verdict: f"{pooled[verdict] / channels:.1%}" for verdict in VERDICTS
}
pd.DataFrame(
    {"channels": pooled, "share": [share_text[verdict] for verdict in VERDICTS]}
).rename(index=plotstyle.VERDICT_LABELS)

WORDS = "zero one two three four five six seven eight nine ten".split()  # for the text

ranked = order.iloc[::-1]  # highest likely share first
top = ranked.iloc[0]
zeros = ranked.loc[ranked["likely_monetized"] == 0, "country"].tolist()
unable_share = snapshot["unable_to_determine"] / checked
unable_share.index = snapshot["country"]
lowest_unable, highest_unable = unable_share.idxmin(), unable_share.idxmax()
print(f"top: {top['country']} {top['likely_share']:.1%}; none likely: {zeros}")
print(f"unable: {unable_share.min():.1%} ({lowest_unable}) "
      f"to {unable_share.max():.1%} ({highest_unable})")

Z = 1.96  # two-sided 95%


def wilson_interval(hits, n):
    """95% Wilson score interval; stays sensible near 0% and at n of about 60."""
    p = hits / n
    centre = (p + Z**2 / (2 * n)) / (1 + Z**2 / n)
    half_width = Z * np.sqrt(p * (1 - p) / n + Z**2 / (4 * n**2)) / (1 + Z**2 / n)
    return centre - half_width, centre + half_width


pooled_likely = pooled["likely_monetized"] / channels
order["low"], order["high"] = wilson_interval(
    order["likely_monetized"], order["channels_checked"]
)
band_low, band_high = wilson_interval(pooled_likely * 60, 60)

outside = order[(order["low"] > pooled_likely) | (order["high"] < pooled_likely)]
expected_outside = 0.05 * countries

# How far apart two countries of 60 channels each must be, near the pooled
# rate, before a two-proportion z-test at the 5% level calls the gap real
needed_gap = Z * np.sqrt(2 * pooled_likely * (1 - pooled_likely) / 60)
neighbour_gaps = ranked["likely_share"].diff().abs().dropna()

print(f"pooled {pooled_likely:.1%}, n=60 band {band_low:.1%} to {band_high:.1%}")
print(f"{len(outside)} of {countries} intervals leave out the pooled rate, "
      f"about {expected_outside:g} expected by chance")
print(f"gap needed {needed_gap:.1%}, "
      f"median gap between neighbours {neighbour_gaps.median():.1%}")

def rank_correlation(a, b):
    """Spearman's rho: Pearson on average ranks, so no SciPy is needed."""
    return float(np.corrcoef(pd.Series(a).rank(), pd.Series(b).rank())[0, 1])


likely, unknown = order["likely_monetized"], order["unable_to_determine"]
order["share_of_called"] = likely / (order["channels_checked"] - unknown)
order["ceiling"] = (likely + unknown) / order["channels_checked"]
order["unable_share"] = unknown / order["channels_checked"]
rho_called = rank_correlation(order["likely_share"], order["share_of_called"])

columns = {
    "country": "country",
    "channels_checked": "checked",
    "unable_share": "unable",
    "likely_share": "likely, of all checked",
    "share_of_called": "likely, of those called",
    "ceiling": "ceiling",
}
(
    order.nlargest(5, "unable_share")[list(columns)]
    .rename(columns=columns)
    .style.format("{:.1%}", subset=list(columns.values())[2:])
    .hide(axis="index")
)

medians = snapshot["median_subscribers"]
over_100k = int((medians > 100_000).sum())
# "all but one country" in the introduction, "49 of the 50" in How a verdict is made
assert over_100k == countries - 1 == 49
rho_size = rank_correlation(order["median_subscribers"], order["likely_share"])
print(f"median channel above 100,000 subscribers in {over_100k} countries; "
      f"median of medians {medians.median():,.0f}; "
      f"rank correlation with likely share {rho_size:.2f}")

bottom_eight = ranked.tail(8)
next_share = ranked.iloc[-9]["likely_share"]
next_up = ranked.loc[ranked["likely_share"] == next_share, "country"].tolist()
print(", ".join(bottom_eight["country"]))
print(f"at most {bottom_eight['likely_share'].max():.1%}; "
      f"next up {next_up} at {next_share:.1%}")