Measles in England and Wales — Analysis¶

# The observed fadeout proportion for each city is the fraction of
# bi-weekly reporting periods with zero reported cases.
obs_prop_zero = (obs_cases == 0).mean(axis=1)   # (954,)

# Group into log10(population) bins for a cleaner look
logpop = np.log10(populations)
print(f"Population range: {populations.min():,} \u2013 {populations.max():,}")
print(f"Observed prop_zero range: {obs_prop_zero.min():.3f} \u2013 {obs_prop_zero.max():.3f}")

nticks = incidence.shape[0]  # 7300
period = 14  # bi-weekly

# Trim to a multiple of 14
n_periods = nticks // period
incidence_trimmed = incidence[:n_periods * period]   # (n_periods*14, 954)
# Reshape and sum: (n_periods, 14, 954) → sum over 14 days → (n_periods, 954)
biweekly_sim = incidence_trimmed.reshape(n_periods, period, -1).sum(axis=1)  # (n_periods, 954)

# Burn in first 2 years (730 days → 52 bi-weekly periods)
burn_periods = (2 * 365) // period
sim_after_burn = biweekly_sim[burn_periods:]   # (n_periods - burn, 954)

sim_prop_zero = (sim_after_burn == 0).mean(axis=0)   # (954,)
print(f"Simulated prop_zero range: {sim_prop_zero.min():.3f} \u2013 {sim_prop_zero.max():.3f}")

# This is the key figure from Conlan et al. (2010):
# smaller cities have higher fadeout fractions.
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(logpop, obs_prop_zero,
           s=18, alpha=0.6, color="#0072B2",
           label="Observed (1944\u20131964)", edgecolors="none")
ax.scatter(logpop, sim_prop_zero,
           s=18, alpha=0.6, color="#E69F00", marker="s",
           label="Simulated (laser.cohorts)", edgecolors="none")

ax.set_xlabel("City population  (log\u2081\u2080 scale)", fontsize=13)
ax.set_ylabel("Proportion of periods\nwith zero reported cases", fontsize=13)
ax.set_title("Measles fadeout and persistence\nEngland & Wales 954 cities, 1944\u20131964", fontsize=13)
ax.set_xticks([3, 4, 5, 6])
ax.set_xticklabels([r"$10^3$", r"$10^4$", r"$10^5$", r"$10^6$"])
ax.legend(fontsize=11)
plt.tight_layout()
plt.show()

highlight = ["London", "Birmingham", "Manchester", "York", "Aberystwyth"]
highlight = [p for p in highlight if p in placenames]

# Observed: bi-weekly reports 1944-1964 (546 periods)
t_obs  = 1944.0 + np.arange(obs_cases.shape[1]) / 26.0  # decimal years

# Simulated: bi-weekly aggregated
n_biweekly_sim = biweekly_sim.shape[0]
t_sim  = 1944.0 + np.arange(n_biweekly_sim) / 26.0

fig, axes = plt.subplots(len(highlight), 2, figsize=(16, 3 * len(highlight)),
                          sharey="row", sharex="col")
fig.suptitle("Observed vs. simulated measles \u2014 selected cities", fontsize=13)

for row, city in enumerate(highlight):
    idx      = placenames.index(city)
    ax_obs   = axes[row, 0]
    ax_sim   = axes[row, 1]

    ax_obs.plot(t_obs, obs_cases[idx], linewidth=0.7, color="#0072B2")
    ax_obs.set_title(f"{city} \u2014 Observed")

    ax_sim.plot(t_sim, biweekly_sim[:, idx], linewidth=0.7, color="#E69F00")
    ax_sim.set_title(f"{city} \u2014 Simulated")

    for ax in (ax_obs, ax_sim):
        ax.set_ylabel("Cases\n(bi-weekly)")

for ax in axes[-1]:
    ax.set_xlabel("Year")

plt.tight_layout()
plt.show()

# Map of England and Wales coloured by fadeout proportion.
# Large cities stay persistent (low prop_zero);
# small peripheral towns fade out often (high prop_zero).
fig, axes = plt.subplots(1, 2, figsize=(14, 10))

for ax, values, title, cmap in [
    (axes[0], obs_prop_zero,  "Observed fadeout proportion",  "Blues"),
    (axes[1], sim_prop_zero, "Simulated fadeout proportion", "Oranges"),
]:
    sc = ax.scatter(longitudes, latitudes,
                    c=values, cmap=cmap, vmin=0, vmax=1,
                    s=4 * np.sqrt(populations / populations.max()) * 80,
                    alpha=0.8, edgecolors="none")
    plt.colorbar(sc, ax=ax, label="Proportion zero")
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.set_title(title)

plt.suptitle("Measles persistence \u2014 England & Wales", fontsize=13)
plt.tight_layout()
plt.show()