nhs-step-up-challenge/Challenge_C_Python_Analysis.py at main · Mulualem03/nhs-step-up-challenge · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
"""
NHS Step Up Career Challenge - Challenge C
Python analysis of longitudinal antidepressant prescribing data, England, Jan 2021 to Oct 2024.

Data source: NHS BSA Open Data Portal, Prescription Cost Analysis (PCA) monthly data
https://opendata.nhsbsa.net/dataset/prescription-cost-analysis-pca-monthly-data

Mulualem, May 2026
"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# ----------------------------------------------------------------------
# 1. LOAD DATA
# ----------------------------------------------------------------------
# Two longitudinal datasets, both from NHS BSA Open Data Portal:
#   - Regional summary: monthly items + cost per region (322 rows)
#   - Regional + drug:  monthly items + cost per region per chemical (9,455 rows)

REG_URL = ('https://raw.githubusercontent.com/nhsengland/Digdata/main/'
           'BSA_ODP_PCA_REGIONAL_SUMMARY.csv')
DRUG_URL = ('https://raw.githubusercontent.com/nhsengland/Digdata/main/'
            'BSA_ODP_PCA_REGIONAL_DRUG_SUMMARY.csv')

reg = pd.read_csv(REG_URL)
drug = pd.read_csv(DRUG_URL)

# Convert YEAR_MONTH (e.g. 202101) to a proper date
reg['date'] = pd.to_datetime(reg['YEAR_MONTH'].astype(str), format='%Y%m')
drug['date'] = pd.to_datetime(drug['YEAR_MONTH'].astype(str), format='%Y%m')

print(f"Regional summary: {reg.shape}, {reg['date'].min():%b %Y} to {reg['date'].max():%b %Y}")
print(f"Drug summary:     {drug.shape}, {drug['BNF_CHEMICAL_SUBSTANCE'].nunique()} drugs")


# ----------------------------------------------------------------------
# 2. NATIONAL-LEVEL HEADLINE NUMBERS
# ----------------------------------------------------------------------
monthly = (reg.groupby('date')
              .agg(items=('ITEMS', 'sum'), cost=('COST', 'sum'))
              .reset_index())

print(f"\nFirst month  (Jan 2021): {monthly['items'].iloc[0]:>12,} items")
print(f"Latest month (Oct 2024): {monthly['items'].iloc[-1]:>12,} items")
growth = (monthly['items'].iloc[-1] / monthly['items'].iloc[0] - 1) * 100
print(f"Growth over the period:    {growth:>+.1f}%")

# Annual totals
reg['year'] = reg['date'].dt.year
yearly = reg.groupby('year').agg(items=('ITEMS', 'sum'),
                                  cost=('COST', 'sum'))
print('\nAnnual totals:')
print(yearly.assign(items_M=yearly['items'] / 1e6,
                    cost_M=yearly['cost'] / 1e6))


# ----------------------------------------------------------------------
# 3. TOP DRUGS NATIONALLY (longitudinal)
# ----------------------------------------------------------------------
top_drugs = (drug.groupby('BNF_CHEMICAL_SUBSTANCE')['ITEMS']
                 .sum()
                 .sort_values(ascending=False)
                 .head(10))
print('\nTop 10 antidepressants over the full period:')
print(top_drugs)


# ----------------------------------------------------------------------
# 4. REGIONAL VARIATION
# ----------------------------------------------------------------------
print('\nRegional growth, Jan 2021 to Oct 2024:')
for region in sorted(reg['REGION_NAME'].unique()):
    if region == 'UNIDENTIFIED':
        continue
    rd = reg[reg['REGION_NAME'] == region].sort_values('date')
    g = (rd['ITEMS'].iloc[-1] / rd['ITEMS'].iloc[0] - 1) * 100
    print(f'  {region:<28s} {rd["ITEMS"].iloc[0]:>10,} -> '
          f'{rd["ITEMS"].iloc[-1]:>10,}  ({g:+.1f}%)')


# ----------------------------------------------------------------------
# 5. FORECAST MODEL
# ----------------------------------------------------------------------
# Linear trend + monthly seasonality dummies. Simple, transparent, and
# easily explainable to non-technical stakeholders.

monthly['t'] = np.arange(len(monthly))
monthly['month'] = monthly['date'].dt.month
month_dummies = pd.get_dummies(monthly['month'], prefix='m',
                                drop_first=True).astype(int)
X = pd.concat([monthly[['t']], month_dummies], axis=1)

m_items = LinearRegression().fit(X, monthly['items'])
m_cost = LinearRegression().fit(X, monthly['cost'])

print(f"\nForecast model fit:")
print(f"  Items R^2 = {m_items.score(X, monthly['items']):.3f}")
print(f"  Cost  R^2 = {m_cost.score(X, monthly['cost']):.3f}")
print(f"  Items trend: {m_items.coef_[0]:+,.0f} items per month")
print(f"  Cost  trend: £{m_cost.coef_[0]:+,.0f} per month")

# Six-month forecast
last_t = monthly['t'].max()
fut_dates = pd.date_range('2024-11-01', periods=6, freq='MS')
fut_t = np.arange(last_t + 1, last_t + 7)
fut_dummies = (pd.get_dummies(fut_dates.month, prefix='m')
                 .reindex(columns=month_dummies.columns, fill_value=0)
                 .astype(int))
X_fut = pd.concat([pd.DataFrame({'t': fut_t}),
                   fut_dummies.reset_index(drop=True)], axis=1)
fc_items = m_items.predict(X_fut)
fc_cost = m_cost.predict(X_fut)

print(f"\nForecast for Nov 2024 to Apr 2025:")
for d, fi, fc in zip(fut_dates, fc_items, fc_cost):
    print(f"  {d:%b %Y}: {fi:>10,.0f} items, £{fc:>13,.0f}")


# ----------------------------------------------------------------------
# 6. VISUALISATIONS
# ----------------------------------------------------------------------
NHS_BLUE = '#005EB8'
NHS_RED = '#DA291C'
plt.rcParams.update({'font.family': 'DejaVu Sans', 'font.size': 11,
                     'axes.spines.top': False, 'axes.spines.right': False,
                     'axes.grid': True, 'grid.alpha': 0.25})

# 6a. National trend with forecast
fig, ax = plt.subplots(figsize=(11, 5))
ax.plot(monthly['date'], monthly['items'] / 1e6,
        color=NHS_BLUE, linewidth=2, marker='o', markersize=3, label='Actual')
ax.plot(fut_dates, fc_items / 1e6,
        color=NHS_RED, linewidth=2, linestyle='--', marker='s',
        markersize=4, label='Forecast')
ax.set_ylabel('Items prescribed (millions)')
ax.set_title('National antidepressant prescribing volume',
             fontweight='bold', loc='left')
ax.legend(frameon=False)
plt.tight_layout()
plt.savefig('national_trend.png', dpi=150, bbox_inches='tight')

# 6b. Volume vs cost (twin axis - the key insight)
fig, ax1 = plt.subplots(figsize=(11, 5))
ax1.plot(monthly['date'], monthly['items'] / 1e6,
         color=NHS_BLUE, linewidth=2, label='Items')
ax1.set_ylabel('Items (millions)', color=NHS_BLUE)
ax1.tick_params(axis='y', labelcolor=NHS_BLUE)

ax2 = ax1.twinx()
ax2.plot(monthly['date'], monthly['cost'] / 1e6,
         color=NHS_RED, linewidth=2, label='Cost')
ax2.set_ylabel('Cost (£ millions)', color=NHS_RED)
ax2.tick_params(axis='y', labelcolor=NHS_RED)
ax2.spines['top'].set_visible(False)
ax2.grid(False)

ax1.set_title('Items rise, cost falls: the generic-drug effect',
              fontweight='bold', loc='left')
plt.tight_layout()
plt.savefig('volume_vs_cost.png', dpi=150, bbox_inches='tight')

print('\nDone. Charts saved.')