Importing Relevant Libraries and dictionaries¶
In [218]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from scipy.stats import chi2_contingency, mannwhitneyu, kruskal, spearmanr
In [219]:
# Display settings
pd.set_option("display.max_columns", None)
pd.set_option("display.max_rows", 100)
pd.set_option("display.width", 120)
pd.set_option("display.float_format", "{:.2f}".format)
# Consistent color palette
COLORS = {
"blue": "#1F77B4",
"orange": "#FF7F0E",
"green": "#2CA02C",
"red": "#D62728",
"purple": "#9467BD",
"grey": "#7F7F7F"
}
# Plot styling
plt.style.use("ggplot")
sns.set_theme(style="whitegrid")
Importing the dataset¶
In [220]:
healthaccess_data = pd.read_csv(r"C:\NYABUTO\DATA ANALYTICS\Data Science and Analytics- Luxdev\Projects\Health Project\Data\hospital_data.csv")
Inspecting the dataset¶
Basic Checks for entire dataset
In [221]:
def data_summary(dataset):
print(f"Dataset summary\n: {dataset.shape}")
print("=" * 50)
print(f"Duplicates \n: {dataset.duplicated().sum()}")
print(f"Columns\n: {dataset.columns}")
print("=" * 50)
print(f"Data Info \n: {dataset.info()}")
print("Checking Missing Values\n", pd.DataFrame({ "Value Count": dataset.isna().sum(), "Percentage": (dataset.isna().sum()/len(dataset)*100).round(2)}))
print("=" * 50)
print(f"Numeric Summary\n: {display(dataset.describe().T)}")
print("=" * 50)
print(f"Categorical Summary: {display(dataset.describe(include="object").T)}")
In [222]:
data_summary(healthaccess_data)
Dataset summary
: (83966, 36)
==================================================
Duplicates
: 0
Columns
: Index(['visit_id', 'first_name', 'last_name', 'gender', 'age', 'county', 'sub_county', 'insurance_type', 'visit_date',
'admission_date', 'discharge_date', 'visit_class', 'visit_type', 'department', 'icd10_code', 'diagnosis_name',
'diagnosis_category', 'triage_level', 'bp_systolic', 'bp_diastolic', 'temperature_c', 'weight_kg', 'spo2_pct',
'length_of_stay_days', 'ward', 'doctor_id', 'referring_facility', 'discharge_outcome', 'readmission_30days',
'bill_amount_ksh', 'sha_covered_ksh', 'payment_method', 'distance_to_facility_km', 'distance_band',
'catchment_zone', 'patient_id'],
dtype='object')
==================================================
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 83966 entries, 0 to 83965
Data columns (total 36 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 visit_id 83966 non-null object
1 first_name 83966 non-null object
2 last_name 83966 non-null object
3 gender 83966 non-null object
4 age 82628 non-null object
5 county 83966 non-null object
6 sub_county 83966 non-null object
7 insurance_type 77931 non-null object
8 visit_date 83966 non-null object
9 admission_date 23983 non-null object
10 discharge_date 23983 non-null object
11 visit_class 83966 non-null object
12 visit_type 83966 non-null object
13 department 83966 non-null object
14 icd10_code 83966 non-null object
15 diagnosis_name 83966 non-null object
16 diagnosis_category 83966 non-null object
17 triage_level 83966 non-null object
18 bp_systolic 83082 non-null object
19 bp_diastolic 83966 non-null int64
20 temperature_c 83966 non-null float64
21 weight_kg 83966 non-null float64
22 spo2_pct 83966 non-null int64
23 length_of_stay_days 83966 non-null int64
24 ward 23983 non-null object
25 doctor_id 83966 non-null object
26 referring_facility 71599 non-null object
27 discharge_outcome 83966 non-null object
28 readmission_30days 83966 non-null object
29 bill_amount_ksh 83966 non-null float64
30 sha_covered_ksh 83966 non-null float64
31 payment_method 83966 non-null object
32 distance_to_facility_km 83966 non-null float64
33 distance_band 83966 non-null object
34 catchment_zone 83966 non-null object
35 patient_id 83966 non-null object
dtypes: float64(5), int64(3), object(28)
memory usage: 23.1+ MB
Data Info
: None
Checking Missing Values
Value Count Percentage
visit_id 0 0.00
first_name 0 0.00
last_name 0 0.00
gender 0 0.00
age 1338 1.59
county 0 0.00
sub_county 0 0.00
insurance_type 6035 7.19
visit_date 0 0.00
admission_date 59983 71.44
discharge_date 59983 71.44
visit_class 0 0.00
visit_type 0 0.00
department 0 0.00
icd10_code 0 0.00
diagnosis_name 0 0.00
diagnosis_category 0 0.00
triage_level 0 0.00
bp_systolic 884 1.05
bp_diastolic 0 0.00
temperature_c 0 0.00
weight_kg 0 0.00
spo2_pct 0 0.00
length_of_stay_days 0 0.00
ward 59983 71.44
doctor_id 0 0.00
referring_facility 12367 14.73
discharge_outcome 0 0.00
readmission_30days 0 0.00
bill_amount_ksh 0 0.00
sha_covered_ksh 0 0.00
payment_method 0 0.00
distance_to_facility_km 0 0.00
distance_band 0 0.00
catchment_zone 0 0.00
patient_id 0 0.00
==================================================
| count | mean | std | min | 25% | 50% | 75% | max | |
|---|---|---|---|---|---|---|---|---|
| bp_diastolic | 83966.00 | 90.01 | 17.66 | 60.00 | 75.00 | 90.00 | 105.00 | 125.00 |
| temperature_c | 83966.00 | 262.54 | 159.41 | 35.50 | 38.70 | 367.00 | 376.00 | 410.20 |
| weight_kg | 83966.00 | 59.32 | 19.43 | 1.20 | 51.50 | 62.70 | 72.10 | 112.10 |
| spo2_pct | 83966.00 | 95.97 | 3.94 | 78.00 | 94.00 | 97.00 | 99.00 | 100.00 |
| length_of_stay_days | 83966.00 | 1.17 | 2.62 | 0.00 | 0.00 | 0.00 | 1.00 | 36.00 |
| bill_amount_ksh | 83966.00 | 15237.00 | 31280.14 | 200.00 | 881.74 | 2376.09 | 15492.60 | 647373.26 |
| sha_covered_ksh | 83966.00 | 2951.29 | 6783.25 | 0.00 | 0.00 | 460.37 | 1461.89 | 83701.16 |
| distance_to_facility_km | 83966.00 | 44.95 | 29.01 | 0.50 | 23.40 | 41.40 | 62.10 | 128.90 |
Numeric Summary : None ==================================================
| count | unique | top | freq | |
|---|---|---|---|---|
| visit_id | 83966 | 83966 | VIS0000001 | 1 |
| first_name | 83966 | 11261 | Akinyi | 488 |
| last_name | 83966 | 102 | Ochieng | 5026 |
| gender | 83966 | 8 | FEMALE | 12876 |
| age | 82628 | 289 | 38.0 | 1829 |
| county | 83966 | 6 | Kisumu | 41667 |
| sub_county | 83966 | 39 | Muhoroni | 6249 |
| insurance_type | 77931 | 7 | No Insurance | 27559 |
| visit_date | 83966 | 1096 | 27/06/2024 | 196 |
| admission_date | 23983 | 1096 | 13/06/2024 | 71 |
| discharge_date | 23983 | 1112 | 26/06/2024 | 74 |
| visit_class | 83966 | 2 | Outpatient | 59983 |
| visit_type | 83966 | 3 | Follow-up | 37547 |
| department | 83966 | 11 | General Practice (GP) | 39518 |
| icd10_code | 83966 | 43 | B20 | 9845 |
| diagnosis_name | 83966 | 43 | HIV/AIDS | 9845 |
| diagnosis_category | 83966 | 16 | Infectious Disease | 17350 |
| triage_level | 83966 | 8 | 5 - Non-Urgent | 22967 |
| bp_systolic | 83082 | 125 | 144 | 875 |
| ward | 23983 | 10 | Resuscitation Bay | 6494 |
| doctor_id | 83966 | 66 | DR065 | 6759 |
| referring_facility | 71599 | 14 | Maseno Mission Hospital | 5535 |
| discharge_outcome | 83966 | 9 | Discharged (Improved) | 34006 |
| readmission_30days | 83966 | 8 | No | 80772 |
| payment_method | 83966 | 10 | SHA | 24946 |
| distance_band | 83966 | 5 | Far (40–80 km) | 33335 |
| catchment_zone | 83966 | 5 | Far | 33335 |
| patient_id | 83966 | 58480 | PAT011414 | 17 |
Categorical Summary: None
Data Cleaning¶
In [223]:
# Drop personally identifiable information (PII) that is not required for analysis.
# Patient names are excluded to protect patient privacy and maintain data confidentiality.
healthaccess_data = healthaccess_data.drop(columns=["first_name","last_name"])
Cleaning Categorical colums¶
In [224]:
categorical_columns = healthaccess_data.select_dtypes(include=["object"]).columns.tolist()
# Removing leading and trailing spaces and case standardization
for column in categorical_columns:
healthaccess_data[column] = healthaccess_data[column].str.strip().str.title()
Gender: Standardizing categories
In [225]:
# Checking for the different variations
healthaccess_data["gender"].unique()
Out[225]:
array(['Female', 'F', 'M', 'Male'], dtype=object)
In [226]:
mapping = {"M": "Male",
"F":"Female"}
healthaccess_data["gender"] = healthaccess_data["gender"].replace(mapping)
Insurance Type
In [227]:
healthaccess_data["insurance_type"].unique()
Out[227]:
array(['Linda Mama', 'No Insurance', 'Sha', 'Corporate',
'Private Insurance', nan], dtype=object)
In [228]:
mapping = {'Sha':'SHA'}
healthaccess_data["insurance_type"] = healthaccess_data["insurance_type"].replace(mapping)
Date Columns
In [229]:
# Identifying entries with both admission dates, discharge dates & Checking date formats
healthaccess_data.loc[healthaccess_data["admission_date"].notnull() & healthaccess_data["discharge_date"].notnull(),["visit_date", "admission_date","discharge_date","length_of_stay_days"]].head(20)
Out[229]:
| visit_date | admission_date | discharge_date | length_of_stay_days | |
|---|---|---|---|---|
| 0 | 21/11/2024 | 21/11/2024 | 25/11/2024 | 4 |
| 20 | 16/12/2023 | 16/12/2023 | 17/12/2023 | 1 |
| 26 | 30/06/2025 | 30/06/2025 | 01/07/2025 | 1 |
| 28 | 17/06/2024 | 17/06/2024 | 22/06/2024 | 5 |
| 34 | 22/08/2025 | 22/08/2025 | 25/08/2025 | 3 |
| 35 | 19/02/2024 | 19/02/2024 | 20/02/2024 | 1 |
| 41 | 19/03/2023 | 19/03/2023 | 20/03/2023 | 1 |
| 42 | 08/11/2024 | 08/11/2024 | 10/11/2024 | 2 |
| 43 | 22/07/2025 | 22/07/2025 | 09/08/2025 | 18 |
| 49 | 09/08/2025 | 09/08/2025 | 11/08/2025 | 2 |
| 52 | 01/09/2024 | 01/09/2024 | 04/09/2024 | 3 |
| 61 | 21/03/2024 | 21/03/2024 | 22/03/2024 | 1 |
| 63 | 02/08/2025 | 02/08/2025 | 03/08/2025 | 1 |
| 66 | 09/10/2024 | 09/10/2024 | 15/10/2024 | 6 |
| 71 | 08/07/2025 | 08/07/2025 | 09/07/2025 | 1 |
| 73 | 03/02/2025 | 03/02/2025 | 04/02/2025 | 1 |
| 76 | 29/11/2023 | 29/11/2023 | 03/12/2023 | 4 |
| 77 | 12/10/2024 | 12/10/2024 | 14/10/2024 | 2 |
| 90 | 12/06/2024 | 12/06/2024 | 18/06/2024 | 6 |
| 91 | 12/01/2025 | 12/01/2025 | 16/01/2025 | 4 |
In [230]:
date_columns = ["visit_date", "admission_date", "discharge_date"]
for dates in date_columns:
healthaccess_data[dates]= pd.to_datetime(healthaccess_data[dates], format = "mixed", errors="coerce")
Triage Level
In [231]:
healthaccess_data["triage_level"].unique()
Out[231]:
array(['4 - Less Urgent', '5 - Non-Urgent', '3 - Urgent', 'Red',
'2 - Emergent', '1 - Resuscitation', 'Yellow', 'Green'],
dtype=object)
In [232]:
mapping = {"Red": "1 - Resuscitation",
"Yellow": "3 - Urgent",
"Green": "5 - Non-Urgent"}
healthaccess_data["triage_level"] = healthaccess_data["triage_level"].replace(mapping)
Referring facility
In [233]:
healthaccess_data["referring_facility"].unique()
Out[233]:
array(['Kisii Teaching & Referral', 'Siaya County Referral',
'Aga Khan Hospital Kisumu', 'Kendu Bay Mission',
'Homabay County Referral', "St. Mary'S Mission Hospital",
'Maseno Mission Hospital', 'Kisumu County Hospital', nan,
'Self Referral', 'Self-Referral', 'Private Clinic',
'Jaramogi Oginga Odinga Th', 'Kakamega County General'],
dtype=object)
In [234]:
mapping = {"Self-Referral": "Self Referral"}
healthaccess_data["referring_facility"] = healthaccess_data["referring_facility"].replace(mapping)
Discharge Outcome
In [235]:
healthaccess_data["discharge_outcome"].unique()
Out[235]:
array(['Discharged (Improved)', 'Discharged (Stable)', 'Transferred',
'Improved', 'Discharged', 'Died', 'Dama', 'Referred Out',
'Discharged - Routine'], dtype=object)
In [236]:
mapping = {"Discharged (Improved)": "Discharged",
"Discharged - Routine": "Discharged",
"Discharged (Stable)": "Discharged",
"Improved": "Discharged",
"Dama": "DAMA"}
healthaccess_data["discharge_outcome"] = healthaccess_data["discharge_outcome"].replace(mapping)
Readmission within 30days
In [237]:
healthaccess_data["readmission_30days"].unique()
Out[237]:
array(['No', 'Yes', 'True', 'False', '1', '0', 'Y'], dtype=object)
In [238]:
mapping = {"Y": "Yes",
"True": "Yes",
"1": "Yes",
"False": "No",
"0":"No"}
healthaccess_data["readmission_30days"] = healthaccess_data["readmission_30days"].replace(mapping)
Payment Method
In [239]:
healthaccess_data["payment_method"].unique()
Out[239]:
array(['Bank Transfer', 'Corporate Cheque', 'Sha', 'Cash', 'Mpesa',
'M-Pesa', 'Linda Mama', 'Waiver/Exemption', 'Sha + Copay'],
dtype=object)
In [240]:
mapping = {"Mpesa": "M-Pesa",
"Sha": "Insurance",
"Sha + Copay": "Insurance + Copay"}
healthaccess_data["payment_method"] = healthaccess_data["payment_method"].replace(mapping)
Cleaning Numeric Columns¶
Age Column
In [241]:
healthaccess_data["age"].unique()
Out[241]:
array(['26.0', '4.0', '34.0', '17.0', '25.0', '55.0', '56.0', '27.0',
'28.0', nan, '18.0', '20.0', '48.0', '29.0', '30.0', '65.0',
'16.0', '22.0', '19.0', '21.0', '58.0', '59.0', '31.0', '32.0',
'33.0', '12.0', '13.0', '14.0', '49.0', '50.0', '38.0', '83.0',
'84.0', '85.0', '1.0', '2.0', '3.0', '7.0', '53.0', '23.0', '57.0',
'0.0', '42.0', '41.0', '43.0', 'Unknown', '9.0', '24.0', '35.0',
'36.0', '60.0', '15.0', '45.0', '39.0', '5.0', '63.0', '76.0',
'22 Yrs', '37.0', '44.0', '80.0', '82.0', '81.0', '46.0', '54.0',
'51.0', '8.0', '52.0', '61.0', '6.0', '42 Yrs', '32 Yrs', '10.0',
'77.0', '78.0', '79.0', '47.0', '75.0', '48 Yrs', '40.0', '69.0',
'70.0', '71.0', '69 Yrs', '11.0', '57 Yrs', '67.0', '68.0',
'51 Yrs', '64.0', '68 Yrs', '72.0', '62.0', '55 Yrs', '49 Yrs',
'66.0', '73.0', '74.0', '66 Yrs', '67 Yrs', '29 Yrs', '43 Yrs',
'10 Yrs', '5 Yrs', '17 Yrs', '44 Yrs', '53 Yrs', '38 Yrs',
'46 Yrs', '34 Yrs', '35 Yrs', '79 Yrs', '27 Yrs', '62 Yrs', '86.0',
'87.0', '15 Yrs', '45 Yrs', '3 Yrs', '31 Yrs', '28 Yrs', '25 Yrs',
'40 Yrs', '24 Yrs', '61 Yrs', '12 Yrs', '23 Yrs', '7 Yrs',
'26 Yrs', '2 Yrs', '11 Yrs', '76 Yrs', '37 Yrs', '36 Yrs',
'59 Yrs', '8 Yrs', '19 Yrs', '21 Yrs', '20 Yrs', '0 Yrs', '39 Yrs',
'14 Yrs', '41 Yrs', '82 Yrs', '54 Yrs', '47 Yrs', '56 Yrs',
'81 Yrs', '16 Yrs', '9 Yrs', '58 Yrs', '4 Yrs', '84 Yrs', '64 Yrs',
'33 Yrs', '18 Yrs', '52 Yrs', '50 Yrs', '63 Yrs', '30 Yrs',
'73 Yrs', '65 Yrs', '70 Yrs', '74 Yrs', '13 Yrs', '80 Yrs',
'60 Yrs', '6 Yrs', '75 Yrs', '78 Yrs', '2.1', '1.9', '3.7', '2.9',
'2.7', '1.1', '1.3', '2.4', '3.2', '0.1', '3.1', '1.8', '4.1',
'1.2', '0.2', '0.7', '4.5', '4.7', '3.3', '0.8', '3.8', '3.4',
'0.5', '2.2', '0.3', '0.9', '3.6', '4.6', '4.8', '1.5', '2.3',
'1.6', '4.4', '4.3', '0.6', '1.4', '1.7', '2.5', '4.2', '0.4',
'2.6', '2.8', '3.9', '4.9', '3.5'], dtype=object)
In [242]:
# Convert unknown to null
healthaccess_data["age"] = healthaccess_data["age"].replace("Unknown", np.nan)
# Remove non numeric entries
healthaccess_data["age"] = healthaccess_data["age"].replace(r" Yrs", "", regex=True)
#Correct datatype
healthaccess_data["age"] = pd.to_numeric(healthaccess_data["age"])
Bp Systolic to numeric dtype
In [243]:
healthaccess_data["bp_systolic"].unique()
Out[243]:
array(['189', '124', '141', '96', '125', '136', '94', '145', '99', '184',
'176', '118', '175', '114', '138', '163', '128', '113', 'Error',
'117', '153', '107', '187', '178', '100', '121', '119', '148',
'126', '139', '161', '135', '116', '179', '115', '122', '188',
'155', nan, '0', '186', '167', '166', '185', '137', '132', '103',
'181', '110', '102', '156', '106', '142', '190', '157', '98', '93',
'182', '120', '158', '91', '177', '133', '109', '95', '127', '147',
'174', '101', '130', '183', '129', '168', '97', '149', '112',
'999', '180', '172', '92', '140', '171', '165', '131', '146',
'164', '105', '143', '108', '104', '151', '162', '169', '159',
'123', '150', '154', '170', '111', '90', '160', '144', '134',
'152', '173', '60', '996', '88', '193', '192', '86', '191', '1004',
'195', '1000', '89', '1003', '997', '87', '85', '998', '995',
'194', '994', '1001', '1002'], dtype=object)
In [244]:
# Convert any text entries to null
healthaccess_data["bp_systolic"] = healthaccess_data["bp_systolic"].replace("Error", np.nan)
#Correct datatype
healthaccess_data["bp_systolic"] = pd.to_numeric(healthaccess_data["bp_systolic"])
Converting Temperature to numeric dtype
In [245]:
# Removing non numeric characters
healthaccess_data["temperature_c"] = healthaccess_data["temperature_c"].replace(r"\D", "",regex=True)
# Correct datatype
healthaccess_data["temperature_c"]= pd.to_numeric(healthaccess_data["temperature_c"])
Handling Missing Values¶
In [246]:
# Share of missing values for columns with missing data
missing_values = (healthaccess_data.isna().sum() / len(healthaccess_data)) * 100
missing_values[missing_values > 0]
Out[246]:
age 4.58 insurance_type 7.19 admission_date 71.44 discharge_date 71.44 bp_systolic 2.06 ward 71.44 referring_facility 14.73 dtype: float64
Age
In [247]:
# Explore if any of the missing are repeat patients that have recorded age in previous visit
missing_patients = healthaccess_data.loc[healthaccess_data["age"].isna(), "patient_id"].unique()
healthaccess_data[healthaccess_data["patient_id"].isin(missing_patients)][["patient_id", "visit_date", "age"]].sort_values(["patient_id", "visit_date"]).head(20)
Out[247]:
| patient_id | visit_date | age | |
|---|---|---|---|
| 10 | Pat000007 | 2023-01-26 | NaN |
| 75292 | Pat000007 | 2023-07-08 | 28.00 |
| 75382 | Pat000007 | 2024-04-20 | 29.00 |
| 8 | Pat000007 | 2024-07-25 | 27.00 |
| 9 | Pat000007 | 2024-11-07 | 28.00 |
| 75569 | Pat000007 | 2025-11-05 | 29.00 |
| 94 | Pat000032 | 2023-08-14 | 43.00 |
| 99 | Pat000032 | 2023-08-23 | NaN |
| 97 | Pat000032 | 2024-11-20 | 43.00 |
| 90 | Pat000032 | 2024-12-06 | 42.00 |
| 98 | Pat000032 | 2025-02-03 | 43.00 |
| 96 | Pat000032 | 2025-08-26 | 43.00 |
| 91 | Pat000032 | 2025-12-01 | 41.00 |
| 155 | Pat000056 | 2023-01-22 | 22.00 |
| 157 | Pat000056 | 2023-05-17 | 22.00 |
| 156 | Pat000056 | 2024-06-26 | NaN |
| 68102 | Pat000065 | 2023-08-08 | 35.00 |
| 190 | Pat000065 | 2024-03-10 | 36.00 |
| 193 | Pat000065 | 2024-12-17 | NaN |
| 189 | Pat000065 | 2025-03-06 | 35.00 |
In [248]:
# Impute with the closest recorded age
healthaccess_data = healthaccess_data.sort_values(["patient_id", "visit_date"])
healthaccess_data["age"] = (healthaccess_data.groupby("patient_id")["age"].ffill().fillna(healthaccess_data.groupby("patient_id")["age"].bfill()))
In [249]:
# Fill the rest with the median age for department
healthaccess_data["age"] = healthaccess_data["age"].fillna(healthaccess_data.groupby("department")["age"].transform("median"))
In [250]:
# Conveting all ages to integers
healthaccess_data["age"] = healthaccess_data["age"].astype(int)
Insurance Type, Triage level, referring facility
In [251]:
missing_categories = ["insurance_type","triage_level", "referring_facility"]
for category in missing_categories:
healthaccess_data[category] = healthaccess_data[category].fillna("Unknown")
Ward
- The same exact number of missing value counts for admission date, dischrge date and ward are expected, because these numbers indicate patients whose visit type was inpatient, who left same day they visited.
In [252]:
# Confirming similar missing entries for admission date, discharge date and ward
healthaccess_data[["admission_date", "discharge_date", "ward"]].isna().value_counts()
Out[252]:
admission_date discharge_date ward True True True 59983 False False False 23983 Name: count, dtype: int64
In [253]:
# Missing values in ward confirmed to be Not applicable(Outpatient visits)
healthaccess_data["ward"]= healthaccess_data["ward"].fillna("Not Applicable")
Data Summary¶
In [254]:
total_visits = len(healthaccess_data)
unique_patients = healthaccess_data['patient_id'].nunique()
total_doctors = healthaccess_data["doctor_id"].nunique()
date_min = healthaccess_data["visit_date"].min()
date_max = healthaccess_data["visit_date"].max()
print("=" * 50)
print("Dataset Summary")
print("=" * 50)
print(f"Total Doctors: {total_doctors:,}")
print(f"Total visits: {total_visits:,}")
print(f"Unique patients: {unique_patients:,}")
print(f"Date range: {date_min.strftime('%d %b %Y')} to {date_max.strftime('%d %b %Y')}")
print()
missing = healthaccess_data.isna().sum()
missing = missing[missing > 0].sort_values(ascending=False)
missing_pct = (missing / total_visits * 100).round(1)
print("Missing data by column:")
print("=" * 50)
for col in missing.index:
print(f" {col:<22} {missing[col]:>6,} missing ({missing_pct[col]}%)")
print()
print("Note: admission_date / discharge_date / ward expected to be blank for outpatient visits")
================================================== Dataset Summary ================================================== Total Doctors: 66 Total visits: 83,966 Unique patients: 58,480 Date range: 01 Jan 2023 to 31 Dec 2025 Missing data by column: ================================================== admission_date 59,983 missing (71.4%) discharge_date 59,983 missing (71.4%) bp_systolic 1,731 missing (2.1%) Note: admission_date / discharge_date / ward expected to be blank for outpatient visits
EDA¶
1. Descriptive Statistics¶
Shape of numeric columns distribution
In [255]:
# Define numeric columns
numerical_columns = healthaccess_data.select_dtypes(include=["int", "float"]).columns.tolist()
In [256]:
# Histograms showing the shape of each distribution
healthaccess_data[numerical_columns].hist(bins = 25, figsize=(10,8))
plt.suptitle("Distribution of each variable")
plt.tight_layout()
plt.show()
Identifying outliers using boxplot
In [257]:
healthaccess_data[numerical_columns].plot(kind='box', subplots=True, layout=(4,4), figsize=(10,8))
plt.suptitle("Distribution of each variable")
plt.tight_layout()
plt.show()
In [258]:
# 1. BP Systolic readings of 0 and 999 are implausible
healthaccess_data.loc[(healthaccess_data["bp_systolic"] == 0) |(healthaccess_data["bp_systolic"] > 250),"bp_systolic"] = np.nan
# 2. Temperature was recorded in tenths of °C; divide by 10 to convert to °C
# First convert datatype to float since it is currently int
healthaccess_data["temperature_c"] = healthaccess_data["temperature_c"].astype(float)
healthaccess_data.loc[healthaccess_data["temperature_c"] > 100, "temperature_c"] = (healthaccess_data.loc[healthaccess_data["temperature_c"] > 100, "temperature_c"] / 10)
In [259]:
# Correct distribution
numerical_columns = healthaccess_data.select_dtypes(include=["int", "float"]).columns.tolist()
healthaccess_data[numerical_columns].hist(bins = 25, figsize=(10,6))
plt.suptitle("Distribution of each variable")
plt.tight_layout()
plt.show()
A. Exploratory Data Analysis (EDA)¶
Patient Demographics
In [260]:
# Patient age distribution
plt.figure(figsize=(4,4))
sns.histplot(data=healthaccess_data["age"].dropna())
plt.title("Patient age distribution")
plt.show()
In [261]:
# Patient gender distribution
gender_counts = healthaccess_data["gender"].value_counts()
labels = []
for gender, count in gender_counts.items():
percentage = count / gender_counts.sum() * 100
labels.append(f"{gender}\n{count:,} ({percentage:.1f}%)")
plt.figure(figsize=(3,4))
plt.pie(gender_counts.values, labels=labels)
plt.title("Share of hospital visits by gender")
plt.show()
In [262]:
# Patient distribution by county (Where are most of the patients coming from?)
county_share = healthaccess_data["county"].value_counts(normalize="percent")*100
county_share
Out[262]:
county Kisumu 49.62 Kakamega 10.89 Migori 10.35 Siaya 10.04 Homa Bay 9.97 Kisii 9.13 Name: proportion, dtype: float64
Patient visits across the years(Hospital_wide capacity)
Annual Visits
In [263]:
yearly_visits = healthaccess_data["visit_date"].dt.year.value_counts()
yearly_visits
Out[263]:
visit_date 2024 31675 2025 30920 2023 21371 Name: count, dtype: int64
In [264]:
# Extracting year and month from visit date column
healthaccess_data["year"] = healthaccess_data["visit_date"].dt.year
healthaccess_data["month"] = healthaccess_data["visit_date"].dt.month
# Visits per year
yearly_visits = healthaccess_data["visit_date"].dt.year.value_counts()
sns.barplot(x = yearly_visits.index, y = yearly_visits.values)
plt.xlabel("Year")
plt.ylabel("Visits")
plt.title("Share of patient visits by year")
plt.show()
Monthly Visits
In [265]:
# Monthly visits over the years
monthly_visits = (healthaccess_data.groupby(["year", "month"]).size().reset_index(name="visits"))
# Visit distribution by year
plt.figure(figsize=(12, 6))
sns.lineplot(data=monthly_visits,x="month",y="visits",hue="year")
plt.xticks(range(1, 13),["Jan", "Feb", "Mar", "Apr", "May", "Jun","Jul", "Aug", "Sep", "Oct", "Nov", "Dec"])
plt.xlabel("Month")
plt.ylabel("Number of Visits")
plt.title("Monthly Hospital Visits by Year")
plt.show()
Daily Visits
In [373]:
# Weekly Visits
healthaccess_data["day_of_week"] = healthaccess_data["visit_date"].dt.day_name()
day_order = ["Monday", "Tuesday", "Wednesday","Thursday", "Friday", "Saturday", "Sunday"]
day_visits = (healthaccess_data["day_of_week"].value_counts().reindex(day_order))
day_visits
Out[373]:
day_of_week Monday 11938 Tuesday 11737 Wednesday 11930 Thursday 12099 Friday 11822 Saturday 12291 Sunday 12149 Name: count, dtype: int64
Number of visits for each visit type
In [267]:
# Visits per month for both inpatient and outpatient
monthly_visit_type = (healthaccess_data.groupby(["year", "month", "visit_class"]).size().reset_index(name="visits"))
# Plot the numbers
monthly_visit_type["date"] = pd.to_datetime(monthly_visit_type[["year", "month"]].assign(day=1))
plt.figure(figsize=(15, 4))
sns.lineplot(data=monthly_visit_type,x="date",y="visits",hue="visit_class")
plt.xlabel("Month")
plt.ylabel("Number of Visits")
plt.title("Monthly Visits by Visit Type")
plt.legend(title="Visit Type")
plt.xticks(monthly_visit_type["date"],monthly_visit_type["date"].dt.strftime("%m/%y"),rotation=45)
plt.tight_layout()
plt.show()
Department - Level Observations
Visits by department
In [268]:
# Number of visits by department
healthaccess_data["department"].value_counts().sort_values(ascending= False)
# Visual
sns.countplot(data =healthaccess_data, x= "department", order=healthaccess_data["department"].value_counts().index)
plt.title("Number of visits by department")
plt.ylabel("Visits")
plt.xlabel("Department")
plt.xticks(rotation = 45, ha = "right")
plt.show()
Inpatient share by department
In [269]:
inpatient_share = healthaccess_data[healthaccess_data["visit_class"]=="Inpatient"]["department"].value_counts(normalize="percent")* 100
# Visual
plt.Figure(figsize=(4,5))
sns.barplot(x=inpatient_share.index, y=inpatient_share.values)
plt.title("Share of inpatient visits by department")
plt.xticks(rotation = 45, ha = "right")
plt.show()
Outpatient share by department
In [270]:
# Department with the highest outpatient visits
outpatient_share = healthaccess_data[healthaccess_data["visit_class"]=="Outpatient"]["department"].value_counts(normalize="percent")*100
# Visual
plt.Figure(figsize=(4,5))
sns.barplot(x= outpatient_share.index, y= outpatient_share.values)
plt.title("Share of outpatient visits by Department")
plt.xticks(rotation = 45, ha = "right")
plt.show()
Top diagnoses in the 2 departments with highest share of visits for inpatient and outpatient(Emergency and casualty for Inpatient and GP for outpatient)
In [271]:
for dept in ["General Practice (Gp)", "Emergency & Casualty"]:
top = healthaccess_data[healthaccess_data['department']==dept]['diagnosis_category'].value_counts(normalize=True).mul(100).round(2). head(2)
print(f"\n{dept}")
print("-" * 60)
for diagnosis, percentage in top.items():
print(f"{diagnosis:<35} {percentage:>6.2f}%")
General Practice (Gp) ------------------------------------------------------------ Parasitic Disease 21.94% Respiratory 15.84% Emergency & Casualty ------------------------------------------------------------ Trauma/Injury 62.01% Parasitic Disease 33.34%
Length of stay by department
In [272]:
# median los
inpatients = healthaccess_data[healthaccess_data["visit_class"] == "Inpatient"]
typical_los_dpt = (inpatients.groupby("department")["length_of_stay_days"].median().sort_values(ascending=False))
typical_los_dpt
Out[272]:
department Icu / Hdu 8.00 Surgery 5.00 Paediatrics 4.00 Oncology 4.00 Tb/Hiv Clinic 3.00 Cardiology 3.00 Maternity 3.00 Emergency & Casualty 2.00 Renal/Dialysis 2.00 Ophthalmology 1.00 Name: length_of_stay_days, dtype: float64
In [273]:
# visual representation
plt.figure(figsize=(10, 6))
sns.boxplot(data=inpatients, x="length_of_stay_days", y="department")
plt.title("Length of Stay by Department")
plt.xlabel("Length of Stay (Days)")
plt.ylabel("Department")
plt.show()
Doctor-patient ratio by department
In [274]:
# Calculating unique patient and doctor IDs by department
unique_patients = healthaccess_data.groupby("department")["patient_id"].nunique()
total_doctors = healthaccess_data.groupby("department")["doctor_id"].nunique()
# Doctor_patient ratio
doctor_patient_ratio = pd.DataFrame({"unique_patients": unique_patients, "total_doctors": total_doctors})
doctor_patient_ratio["patients_per_doctor"] = (doctor_patient_ratio["unique_patients"] /doctor_patient_ratio["total_doctors"]).round()
doctor_patient_ratio = doctor_patient_ratio.sort_values("patients_per_doctor",ascending=False)
doctor_patient_ratio
Out[274]:
| unique_patients | total_doctors | patients_per_doctor | |
|---|---|---|---|
| department | |||
| General Practice (Gp) | 32972 | 11 | 2997.00 |
| Tb/Hiv Clinic | 8821 | 7 | 1260.00 |
| Maternity | 7019 | 7 | 1003.00 |
| Emergency & Casualty | 6808 | 9 | 756.00 |
| Oncology | 3391 | 5 | 678.00 |
| Renal/Dialysis | 3091 | 6 | 515.00 |
| Cardiology | 2856 | 6 | 476.00 |
| Icu / Hdu | 2344 | 7 | 335.00 |
| Ophthalmology | 995 | 3 | 332.00 |
| Surgery | 2173 | 7 | 310.00 |
| Paediatrics | 214 | 5 | 43.00 |
Referrals received by department
In [275]:
# Share of referrals by dpt
referral_by_dpt = healthaccess_data[healthaccess_data["visit_type"]=="Referral In"]["department"].value_counts(normalize = "percent")*100
referral_by_dpt
Out[275]:
department Emergency & Casualty 19.83 Tb/Hiv Clinic 17.33 Renal/Dialysis 11.91 General Practice (Gp) 10.83 Oncology 10.19 Maternity 7.79 Icu / Hdu 7.29 Cardiology 6.74 Surgery 6.60 Ophthalmology 1.15 Paediatrics 0.35 Name: proportion, dtype: float64
Readmission rate by department
In [276]:
crosstab = (pd.crosstab(healthaccess_data["department"],healthaccess_data["readmission_30days"], normalize = "index")*100)
crosstab
Out[276]:
| readmission_30days | No | Yes |
|---|---|---|
| department | ||
| Cardiology | 95.50 | 4.50 |
| Emergency & Casualty | 93.34 | 6.66 |
| General Practice (Gp) | 98.21 | 1.79 |
| Icu / Hdu | 93.88 | 6.12 |
| Maternity | 96.52 | 3.48 |
| Oncology | 95.21 | 4.79 |
| Ophthalmology | 98.90 | 1.10 |
| Paediatrics | 96.40 | 3.60 |
| Renal/Dialysis | 98.62 | 1.38 |
| Surgery | 94.48 | 5.52 |
| Tb/Hiv Clinic | 96.22 | 3.78 |
Patient Outcomes
In [277]:
# New Vs Repeat patients
patient_visits = healthaccess_data["patient_id"].value_counts()
repeat_patients = (patient_visits > 1).sum()
unique_patients = patient_visits.count()
print(repeat_patients, unique_patients)
8981 58480
Repeat rates by diagnosis
In [278]:
# No of times each patient was seen for the SAME diagnosis
diag_visit_counts = healthaccess_data.groupby(["patient_id", "diagnosis_name"]).size()
repeat_patients = diag_visit_counts[diag_visit_counts > 1]
repeat_rate = (repeat_patients.groupby("diagnosis_name").size() / diag_visit_counts.groupby("diagnosis_name").size()* 100).sort_values(ascending=False).head(10)
ax = sns.barplot(x=repeat_rate.values, y=repeat_rate.index)
plt.title("Top 10 Diagnoses by Repeat-Patient Rate (%)")
plt.xlabel("Share of Patients with Repeat Visits (%)")
plt.show()
Admission rates
In [279]:
# Proportion of visits resulting in admission
admission_percentage = (((healthaccess_data["visit_class"] == "Inpatient").sum()/ len(healthaccess_data)) * 100).round()
admission_percentage
Out[279]:
np.float64(29.0)
In [280]:
# To find out the share of visits that end up as admitted or outpatient
pd.crosstab(healthaccess_data["visit_type"],healthaccess_data["visit_class"],normalize="index") * 100
Out[280]:
| visit_class | Inpatient | Outpatient |
|---|---|---|
| visit_type | ||
| Follow-Up | 18.46 | 81.54 |
| New/Emergency | 29.94 | 70.06 |
| Referral In | 57.84 | 42.16 |
Discharge Outcomes
In [281]:
# Patient outcomes
healthaccess_data["discharge_outcome"].value_counts(normalize="percent")*100
plt.figure(figsize=(4,4))
sns.countplot(data = healthaccess_data, x = "discharge_outcome", order = healthaccess_data["discharge_outcome"].value_counts().index)
plt.title("Distribution of patient discharge types")
plt.xticks(rotation = 45, ha = "right")
plt.show()
In [282]:
# Discharge outcome by triage level
outcome_by_triage = pd.crosstab(healthaccess_data['triage_level'], healthaccess_data['discharge_outcome'], normalize='index').mul(100).round(1)
print(outcome_by_triage)
discharge_outcome DAMA Died Discharged Referred Out Transferred triage_level 1 - Resuscitation 19.00 15.50 43.00 8.70 13.80 2 - Emergent 10.30 7.50 58.90 10.20 13.10 3 - Urgent 5.30 2.20 74.40 9.90 8.20 4 - Less Urgent 2.20 1.20 85.20 7.20 4.20 5 - Non-Urgent 1.40 0.90 88.60 6.80 2.40
In [283]:
# Patient readmission after 30 days
crosstab = pd.crosstab(healthaccess_data["visit_class"],healthaccess_data["readmission_30days"], normalize = "index")*100
crosstab
Out[283]:
| readmission_30days | No | Yes |
|---|---|---|
| visit_class | ||
| Inpatient | 93.21 | 6.79 |
| Outpatient | 98.35 | 1.65 |
Service Payment & Financial Performance
Insurance type distribution
In [284]:
#Insurance/ payment type distribution
insurance_type = healthaccess_data[healthaccess_data["insurance_type"] != "Unknown"]["insurance_type"].value_counts(normalize="percent")*100
insurance_type
Out[284]:
insurance_type SHA 41.89 No Insurance 35.36 Private Insurance 14.34 Corporate 6.13 Linda Mama 2.28 Name: proportion, dtype: float64
Bill amount distribution
In [285]:
# Bill amount distribution
sns.kdeplot(data=healthaccess_data.dropna(subset=["bill_amount_ksh"]), x="bill_amount_ksh")
plt.title("Distribution of Hospital Bill Amounts")
plt.xlabel("Bill Amount (KSH)")
plt.ylabel("Density")
plt.show()
print("Typical bill amount : Kshs.", healthaccess_data["bill_amount_ksh"].median().round(2))
Typical bill amount : Kshs. 2376.09
Average bill amount for inpatient & Outpatient
In [286]:
avg_bill_visit_class = healthaccess_data.groupby("visit_class")["bill_amount_ksh"].median().sort_values(ascending = False)
avg_bill_visit_class
Out[286]:
visit_class Inpatient 28948.69 Outpatient 1184.94 Name: bill_amount_ksh, dtype: float64
Comparison between the total bill by dpt and average bill(Departments with highest revenues for hospital and highest financial burden to patients)
In [287]:
# Function that relabels axes appropriately
def relabel_axis(ax):
ticks = ax.get_xticks()
ax.set_xticks(ticks)
labels = []
for t in ticks:
if t >= 1_000_000:
labels.append(f'{t/1_000_000:.1f}M')
elif t >= 10_000:
labels.append(f'{t/1_000:.0f}K')
elif t >= 1_000:
labels.append(f'{t/1_000:.1f}K')
else:
labels.append(f'{t:.0f}')
ax.set_xticklabels(labels)
# Total bill by department
bill_amount_dpt = healthaccess_data.groupby("department")["bill_amount_ksh"].sum().sort_values(ascending=False)
ax = sns.barplot(x=bill_amount_dpt.values, y=bill_amount_dpt.index)
relabel_axis(ax)
plt.title("Total Bill Amount by Department")
plt.show()
# Typical bill by department
avg_bill_dpt = healthaccess_data.groupby("department")["bill_amount_ksh"].median().sort_values(ascending=False)
ax = sns.barplot(x=avg_bill_dpt.values, y=avg_bill_dpt.index)
relabel_axis(ax)
plt.title("Average Bill Amount by Department")
plt.show()
# Share of total bill covered by SHA, by department
sha_paid_dpt = healthaccess_data.groupby("department")["sha_covered_ksh"].sum()
sha_share_dpt = (sha_paid_dpt / bill_amount_dpt * 100).sort_values(ascending=False)
ax = sns.barplot(x=sha_share_dpt.values, y=sha_share_dpt.index, order=sha_share_dpt.index)
plt.title("Share of Total Bill Covered by SHA by Department")
plt.xlabel("SHA Coverage (%)")
plt.show()
In [288]:
# SHA coverage by Insurance types
# Create a coverage column
healthaccess_data["sha_coverage_pct"] = (healthaccess_data["sha_covered_ksh"] /healthaccess_data["bill_amount_ksh"]) * 100
coverage_by_insurance = ( healthaccess_data[healthaccess_data["insurance_type"] != "Unknown"].groupby("insurance_type")["sha_coverage_pct"].median().sort_values(ascending=False))
coverage_by_insurance
Out[288]:
insurance_type SHA 73.08 Linda Mama 45.45 Private Insurance 16.04 Corporate 12.39 No Insurance 0.00 Name: sha_coverage_pct, dtype: float64
2. Inferrential Statistics¶
Correlation Analysis¶
In [289]:
correlation = healthaccess_data[numerical_columns].corr(method = "spearman")
plt.figure(figsize = (10,5))
sns.heatmap(data = correlation , annot = True, fmt = ".2f")
plt.show()
Inspecting correlation between numeric columns with relationship for effect size and p values
Further inspecting SHA covered and Length of Stay(I suspect it is a result of the bill amount)
In [290]:
# Spearman correlations and p-values
r_los_sha, p_los_sha = spearmanr(healthaccess_data["length_of_stay_days"],healthaccess_data["sha_covered_ksh"])
r_los_bill, p_los_bill = spearmanr(healthaccess_data["length_of_stay_days"],healthaccess_data["bill_amount_ksh"])
r_sha_bill, p_sha_bill = spearmanr(healthaccess_data["sha_covered_ksh"],healthaccess_data["bill_amount_ksh"])
r_age, p_age = spearmanr(healthaccess_data["age"],healthaccess_data["weight_kg"])
# Partial correlation
partial_r = (r_los_sha - (r_los_bill * r_sha_bill)) / (
((1 - r_los_bill**2) ** 0.5) *((1 - r_sha_bill**2) ** 0.5))
# Coefficients of determination
r2_los_sha = r_los_sha**2
r2_los_bill = r_los_bill**2
r2_sha_bill = r_sha_bill**2
r2_age = r_age**2
r2_partial = partial_r**2
print(f"LOS vs SHA coverage: r = {r_los_sha:.3f}, p = {p_los_sha:.4f}, R² = {r2_los_sha:.3f}")
print()
print(f"LOS vs Bill amount: r = {r_los_bill:.3f}, p = {p_los_bill:.4f}, R² = {r2_los_bill:.3f}")
print()
print(f"SHA coverage vs Bill amount: r = {r_sha_bill:.3f}, p = {p_sha_bill:.4f}, R² = {r2_sha_bill:.3f}")
print()
print(f"Age vs Weight: r = {r_age:.3f}, p = {p_age:.4f}, R² = {r2_age:.3f}")
print()
print(f"Partial LOS vs SHA coverage: r = {partial_r:.3f}, R² = {r2_partial:.3f}")
LOS vs SHA coverage: r = 0.305, p = 0.0000, R² = 0.093 LOS vs Bill amount: r = 0.698, p = 0.0000, R² = 0.487 SHA coverage vs Bill amount: r = 0.420, p = 0.0000, R² = 0.176 Age vs Weight: r = 0.299, p = 0.0000, R² = 0.090 Partial LOS vs SHA coverage: r = 0.018, R² = 0.000
Strong relationship: Length of stay ↔ bill amount - more days admitted, more charges. Weaker positives are mostly downstream of this: bill ↔ SHA coverage (coverage is calculated from bill size), SHA coverage ↔ length of stay (likely mediated through bill amount, not direct), and weight ↔ age (expected physiological pattern).
Hypothesis Testing¶
1. Categorical vs Categorical Variables¶
Given the size of the dataset, the effect size will accompany the p value to determine strength of association
In [291]:
#Cramers V function
def cramers_v(chi2, table):
n = table.sum().sum()
r, c = table.shape
return np.sqrt(chi2 / (n * min(r - 1, c - 1)))
i). Do some departments show significantly higher death rates than others?
In [292]:
# Column with patients outcomes (Dead and not dead)
healthaccess_data["patient_outcome"] = np.where(healthaccess_data["discharge_outcome"] == "Died","Died","Not Died")
In [293]:
# # Chi2 test for the patient outcome vs department
H0 = "There is no association between department and death outcome."
H1 = "The proportion of patients who die differs across departments."
table = pd.crosstab(healthaccess_data["department"], healthaccess_data["patient_outcome"])
α = 0.05
chi2, p, dof, expected = chi2_contingency(table)
print(f"Test statistic (Chi2): {chi2:.2f}")
print(f"P_value : {p:.3f}")
print(f"Degrees of freedom: {dof}")
if p <= α:
print(f"Reject the null hypothesis: {H0}")
else:
print(f"Fail to reject the null hypothesis {H0}")
Cramer_V = cramers_v(chi2, table)
print(f"Effect Size: {Cramer_V:.2f}")
Test statistic (Chi2): 1372.07 P_value : 0.000 Degrees of freedom: 10 Reject the null hypothesis: There is no association between department and death outcome. Effect Size: 0.13
- From the Effect Size 0.13, we can deduce that although the result is statistically significant, the death outcomes and departments are very weakly associated, showing it lacks practical importance.
ii). Is mortality higher among patients classified as more urgent?
In [294]:
table = pd.crosstab(healthaccess_data["triage_level"],healthaccess_data["patient_outcome"])
table
H0 = "Mortality rates are the same among triage levels"
H1 = "Mortality rates differ significantly among triage levels"
α = 0.05
chi2, p, dof, expected = chi2_contingency(table)
print(f"Chi2: {chi2:.3f}")
print(f"Degrees of Freedom: {dof}")
print(f"P_Value: {p:2f}")
if p <= α:
print("Reject the null hypothesis")
else:
print("Fail to reject null Hypothesis")
Cramer_V = cramers_v(chi2, table)
print(f"Effect Size: {Cramer_V:.2f}")
Chi2: 4940.033 Degrees of Freedom: 4 P_Value: 0.000000 Reject the null hypothesis Effect Size: 0.24
iii). Does a patient's diagnosis category influence whether they are admitted as an inpatient or managed as an outpatient?
In [295]:
table = pd.crosstab(healthaccess_data["diagnosis_category"], healthaccess_data["visit_class"])
H0 = "Patient diagnosis has no effect on their admission status."
H1 = "Patient diagnosis significantly affects admission status"
α = 0.05
chi2, p, dof, expected = chi2_contingency(table)
print(f"chi2_contingency: {chi2:.3f}")
print(f"Degrees of Freedom: {dof}")
print(f"P_Value: {p:2f}")
if p <= α:
print(f"Reject the null hypothesis: {H0}")
else:
print(f"Fail to reject null Hypothesis: {H0}")
Cramer_V = cramers_v(chi2, table)
print(f"Effect Size: {Cramer_V:.2f}")
chi2_contingency: 30150.525 Degrees of Freedom: 15 P_Value: 0.000000 Reject the null hypothesis: Patient diagnosis has no effect on their admission status. Effect Size: 0.60
- The effect size of 0.6 means that the patient diagnosis and admission status are strongly associated
iv). Is there a relationship between a patient's gender and the hospital department they visit?
Maternity retains a small share of male patients (~3.5%) from child growth monitoring visits, where the infant rather than the mother is the recorded patient, this doesn't change the department's overall female-skew, so it is excluded from the test.
In [296]:
H0 = "Patient gender does not determine the department they visit."
H1 = "Patient gender determines the department they visit."
table1 = pd.crosstab(healthaccess_data["gender"], healthaccess_data["department"])
table2 = table1.drop(columns = "Maternity")
α = 0.05
chi2, p, dof, expected = chi2_contingency(table1)
chi2_2, p_2, dof_2, expected_2 = chi2_contingency(table2)
print(f"chi2_with Maternity: {chi2:<22.3f} | chi2_without Maternity: {chi2_2:.3f}")
print(f"P_Value with Maternity: {p:<19} | P_Value without Maternity: {p_2}")
print(f" Effect Size With Maternity: {cramers_v(chi2, table1):>14.2f} | Effect Size (Without Maternity):{cramers_v(chi2_2, table2):.2f}")
if p <= α and p_2 <= α:
print()
print(f"Reject the null hypothesis")
else:
print(f"Fail to reject null Hypothesis")
print()
chi2_with Maternity: 4585.261 | chi2_without Maternity: 336.279 P_Value with Maternity: 0.0 | P_Value without Maternity: 5.141883368630628e-67 Effect Size With Maternity: 0.23 | Effect Size (Without Maternity):0.07 Reject the null hypothesis
The gender - department association was small-to-moderate overall (V = 0.23), but dropped to negligible (V = 0.07) once Maternity was excluded. This shows the association is driven almost entirely by Maternity's female skew, not a broader gender pattern across the hospital.
v). Are admitted patients more likely to pay out-of-pocket vs. insurance?
In [297]:
payment_groups = {
"Bank Transfer": "Out-of-pocket",
"M-Pesa": "Out-of-pocket",
"Corporate Cheque": "Out-of-pocket",
"Cash": "Out-of-pocket",
"Insurance": "Insurance",
"Insurance + Copay": "Insurance",
"Linda Mama": "Insurance",
"Waiver/Exemption": "Waiver/Exemption"}
healthaccess_data["payment_group"] = healthaccess_data["payment_method"].map(payment_groups)
In [298]:
H0 = "Visit class does not determine payment group ."
H1 = "Visit class determines payment group."
table1 = pd.crosstab(healthaccess_data["visit_class"], healthaccess_data["payment_group"])
table2 = table1.drop(columns = "Waiver/Exemption")
α = 0.05
chi2, p, dof, expected = chi2_contingency(table1)
chi2_2, p_2, dof_2, expected_2 = chi2_contingency(table2)
print(f"chi2_with Waiver: {chi2:<25.3f} | chi2_without Waiver: {chi2_2:.3f}")
print()
print(f"P_Value with Waiver: {p:<22} | P_Value without Waiver: {p_2}")
print()
print(f"Effect Size (With Waiver): {cramers_v(chi2, table1):<16.2f} | Effect Size (Without Waiver): {cramers_v(chi2_2, table2):.2f}")
print()
if p <= α and p_2 <= α:
print(f"Reject the null hypothesis: {H0}")
else:
print(f"Fail to reject null Hypothesis: {H0}")
chi2_with Waiver: 1.751 | chi2_without Waiver: 0.053 P_Value with Waiver: 0.41673251141121614 | P_Value without Waiver: 0.81720106261047 Effect Size (With Waiver): 0.00 | Effect Size (Without Waiver): 0.00 Fail to reject null Hypothesis: Visit class does not determine payment group .
- No significant association was found between visit class and payment group (with Waiver: χ² = 1.75, p = 0.42; without: χ² = 0.05, p = 0.82), and effect sizes were negligible (V = 0.00) in both cases. Admission status does not predict how patients pay.
2. Categorical vs Numeric variables¶
i). Does a patient's age determine the mortality probability?
Given the size of the dataset, the effect size will accompany the p value to determine strength of association
In [299]:
def rank_biserial(stat, n1, n2):
return 1 - (2 * stat) / (n1 * n2)
In [300]:
group_died = healthaccess_data.loc[healthaccess_data["patient_outcome"] == "Died", "age"]
group_not_died = healthaccess_data.loc[healthaccess_data["patient_outcome"] == "Not Died", "age"]
H0 = "There is no difference in age between patients who died and those who did not."
H1 = "There is a difference in age between patients who died and those who did not."
α = 0.05
stat, p = mannwhitneyu(group_died, group_not_died)
print(f"Mann-Whitney U Statistic: {stat:.3f}")
print(f"P_Value: {p:2f}")
if p <= α:
print(f"Reject the null hypothesis: {H0}")
else:
print(f"Fail to reject null Hypothesis: {H0}")
n1 = len(group_died)
n2 = len(group_not_died)
effect_size = rank_biserial(stat, n1, n2)
print(f"Effect Size (rank-biserial r): {effect_size:.2f}")
print("Median age (Died):", group_died.median())
print("Median age (Not Died):", group_not_died.median())
Mann-Whitney U Statistic: 121312203.000 P_Value: 0.000260 Reject the null hypothesis: There is no difference in age between patients who died and those who did not. Effect Size (rank-biserial r): -0.04 Median age (Died): 37.0 Median age (Not Died): 36.0
- While the mannwhitney U test found a statistically significant association between age and the patient outcomes, the effect size (-0.04) indicates a negligible association, suggesting that the difference is statistically significant but very small in practical terms.
ii). Does the bill amount differ between insured and uninsured patients?
In [301]:
# Create a category that has insured and uninsured
healthaccess_data["Insurance_Status"] = np.where(healthaccess_data["insurance_type"]=="No Insurance","Uninsured", "Insured")
In [302]:
group_insured = healthaccess_data.loc[healthaccess_data["Insurance_Status"]=="Insured", "bill_amount_ksh"]
group_uninsured = healthaccess_data.loc[healthaccess_data["Insurance_Status"]=="Uninsured", "bill_amount_ksh"]
H0 = "The typical bill amount for insured and uninsured patients is the same"
H1 = "The typical bill amount differs between insured and uninsured patients"
α = 0.05
stat, p = mannwhitneyu(group_insured, group_uninsured, alternative="two-sided")
print(f"Mann-Whitney U Statistic: {stat:.3f}")
print(f"P_Value: {p:2f}")
if p <= α:
print("Reject the null Hypothesis")
else:
print("Fail to reject the null hypothesis")
n1 = len(group_insured)
n2 = len(group_uninsured)
effect_size = rank_biserial(stat, n1, n2)
print(f"Effect Size (rank-biserial r): {effect_size:.2f}")
print("Median bill (Insured):", group_insured.median())
print("Median bill (Not insured):", group_uninsured.median())
Mann-Whitney U Statistic: 792552287.500 P_Value: 0.000004 Reject the null Hypothesis Effect Size (rank-biserial r): -0.02 Median bill (Insured): 2476.24 Median bill (Not insured): 2199.08
- The Mann–Whitney U test found a statistically significant difference in bill amounts between insured and uninsured patients. Insured patients had a higher median bill (KSh 2,476.24) than uninsured patients (KSh 2,199.08). However, the negligible effect size (r = −0.02) indicates a negligible effect, suggesting the difference is very small in practical terms.
iii). Does the distance to facility influence patient outcomes?
In [303]:
group_died = healthaccess_data.loc[healthaccess_data["patient_outcome"] == "Died", "distance_to_facility_km"]
group_not_died = healthaccess_data.loc[healthaccess_data["patient_outcome"] == "Not Died", "distance_to_facility_km"]
H0 = "The distance to the health facility does not affect the patient outcomes."
H1 = "The distance to the health facility significantly affects patient outcomes."
α = 0.05
stat, p = mannwhitneyu(group_died, group_not_died)
print(f"Mann-Whitney U Statistic: {stat:.3f}")
print(f"P_Value: {p:2f}")
if p <= α:
print(f"Reject the null hypothesis: The distance has an effect on patient outcomes")
else:
print(f"Fail to reject null Hypothesis: The distance has no effect on patient outcomes")
n1 = len(group_died)
n2 = len(group_not_died)
effect_size = rank_biserial(stat, n1, n2)
print(f"Effect Size (rank-biserial r): {effect_size:.2f}")
print("Median distance (Died):", group_died.median())
print("Median distance (Not Died):", group_not_died.median())
Mann-Whitney U Statistic: 115338311.500 P_Value: 0.305852 Fail to reject null Hypothesis: The distance has no effect on patient outcomes Effect Size (rank-biserial r): 0.01 Median distance (Died): 40.9 Median distance (Not Died): 41.4
The Mann Whitney U test found no significant difference in distance to the facility between patients who died and those who did not (p = 0.306). Median distances were similar (40.9 km vs. 41.4 km), with a negligible effect size (r = 0.01).
iv). Does length of stay differ by department?
In [304]:
def epsilon_squared(stat, n, k):
return (stat - k + 1) / (n - k)
In [305]:
H0 = "The average length of days stayed is the same across departments"
H1 = "The average length of stay differs across departments"
α = 0.05
groups= []
for dpt in healthaccess_data["department"].unique():
days = healthaccess_data[healthaccess_data["department"] == dpt]["length_of_stay_days"]
groups.append(days)
stat, p = kruskal(*groups)
print(f"Kruskal-Wallis H: {stat:.2f}")
print(f"p-value: {p:.4f}")
if p <= α:
print("Reject the null hypothesis: The average LOS differs across departments")
else:
print("Fail to reject the null hypothesis: The average LOS across departments is the same")
n = len(healthaccess_data)
k = healthaccess_data["department"].nunique()
effect_size = epsilon_squared(stat, n, k)
print(f"Effect Size (epsilon-squared): {effect_size:.3f}")
Kruskal-Wallis H: 52188.22 p-value: 0.0000 Reject the null hypothesis: The average LOS differs across departments Effect Size (epsilon-squared): 0.622
- Length of stay differed significantly across departments (H = 52188.22, p < .001, ε² = 0.62) - a large effect, though expected since several departments are inherently inpatient (Surgery, ICU) or outpatient (GP, Ophthalmology) by nature.
v). Triage level vs bill amount
In [306]:
H0 = "The distribution of bill amounts is the same across triage levels."
H1 = "The distribution of bill amounts differs across triage levels."
α = 0.05
groups= []
for level in healthaccess_data["triage_level"].unique():
bill = healthaccess_data[healthaccess_data["triage_level"] == level]["bill_amount_ksh"]
groups.append(bill)
stat, p = kruskal(*groups)
print(f"Kruskal-Wallis H: {stat:.2f}")
print(f"p-value: {p:.4f}")
if p <= α:
print("Reject the null hypothesis: Bill amounts differ across triage levels ")
else:
print("Fail to reject the null hypothesis: Bill amounts are the same across triage levels")
n = len(healthaccess_data)
k = healthaccess_data["triage_level"].nunique()
effect_size = epsilon_squared(stat, n, k)
print(f"Effect Size (epsilon-squared): {effect_size:.3f}")
Kruskal-Wallis H: 1392.42 p-value: 0.0000 Reject the null hypothesis: Bill amounts differ across triage levels Effect Size (epsilon-squared): 0.017
- There is a statistically significant difference in bill amounts across triage levels (H = 1392.42, p < 0.001), but the effect size is small (ε² = 0.017), indicating that triage level explains only a small proportion of the variation in bill amounts.