Psychological Distress and Physical Multimorbidity in England: Evidence from the Health Survey for England 2022¶
Healthcare and Medical Analytics – Group Project¶
Dataset: Health Survey for England 2022 (HSE 2022), UK Data Service, Study Number 9469. NHS England / UKDS.¶
Note on Code Structure: The sections in this Jupyter Notebook are labelled to correspond to specific appendix sections in the main report (e.g. Appendix A.1, Appendix B.2). These labels are used for navigation purposes only and indicate where the corresponding narrative and output appear in the report appendix. The code sections do not represent the appendix content itself, but provide the underlying analysis that supports it.
Note on AI Usage: This analysis was conducted with the assistance of AI tools. Claude (Anthropic) and ChatGPT (OpenAI) were used to support code generation, debugging, and analytical decisions throughout this project. We remain responsible for all analytical choices, interpretations, and conclusions presented in the report. All AI-generated code was reviewed, tested, and validated by myself prior to inclusion.
Imports, Loading, and Cleaning¶
In [5]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
import seaborn as sns
import statsmodels.api as sm
from statsmodels.discrete.discrete_model import NegativeBinomial
from scipy.stats import chi2
import warnings
warnings.filterwarnings('ignore')
In [7]:
df_raw = pd.read_stata('Data/stata/stata13/hse_2022_eul_v1.dta', convert_categoricals=False)
print(f"Dataset loaded successfully.")
print(f"Observations: {df_raw.shape[0]}")
print(f"Variables: {df_raw.shape[1]}")
Dataset loaded successfully. Observations: 9122 Variables: 1883
In [8]:
def clean(series, extra_missing=None):
"""
Replace HSE missing codes (all negative values) with NaN.
HSE uses -1 (not applicable), -8 (don't know), -9 (refused).
extra_missing: additional values to treat as missing
"""
s = series.copy().astype(float)
s[s < 0] = np.nan
if extra_missing:
for v in extra_missing:
s[s == v] = np.nan
return s
df = df_raw.copy()
# Health measures
df['ghq12scr'] = clean(df['ghq12scr'])
df['ghqg2'] = clean(df['ghqg2'])
df['ghq'] = clean(df['ghq'])
df['condlcnt'] = clean(df['condlcnt'])
df['condlcnt2'] = clean(df['condlcnt2'])
df['limlast'] = clean(df['limlast'])
# Key moderator
df['qimd19'] = clean(df['qimd19'])
# Demographics
df['Sex'] = clean(df['Sex'])
df['ag16g10'] = clean(df['ag16g10'])
df['origin2'] = clean(df['origin2'])
# Socioeconomic
df['eqv5'] = clean(df['eqv5'], extra_missing=[-90])
df['nssec3'] = clean(df['nssec3'], extra_missing=[99])
df['topqual3'] = clean(df['topqual3'])
# Lifestyle
df['cigst1_19'] = clean(df['cigst1_19'])
df['totalwug2_22'] = clean(df['totalwug2_22'])
# Geography
df['GOR1'] = clean(df['GOR1'])
df['urban14b'] = clean(df['urban14b'])
# Survey weights
df['wt_int'] = clean(df['wt_int'])
# Physical condition flags (descriptives only)
for v in ['compexp1','compexp3','compexp4',
'compexp6','compexp7','compexp9']:
df[v] = clean(df[v])
# BMI (Appendix A.8)
df['bmisrg3'] = clean(df['bmisrg3'])
print("Missing value coding complete.")
Missing value coding complete.
In [9]:
n_before = len(df)
df = df[df['ag16g10'].notna()].copy()
n_after = len(df)
print(f"Full sample: {n_before:,}")
print(f"Adults 16+ only: {n_after:,}")
print(f"Children excluded: {n_before - n_after:,}")
Full sample: 9,122 Adults 16+ only: 7,729 Children excluded: 1,393
In [50]:
# ------------------------------------------------
# Global plot style
# ------------------------------------------------
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
import numpy as np
import pandas as pd
plt.rcParams.update({
"figure.dpi": 150,
"savefig.dpi": 300,
"font.size": 11,
"axes.titlesize": 13,
"axes.labelsize": 11,
"legend.fontsize": 10,
"legend.title_fontsize": 10,
"xtick.labelsize": 10,
"ytick.labelsize": 10,
"axes.spines.top": False,
"axes.spines.right": False
})
ghq_palette = {
"GHQ 0 (No distress)": "#2C7FB8",
"GHQ 1-3 (Mild)": "#7FCDBB",
"GHQ 4+ (Probable distress)": "#D95F0E"
}
ghq_markers = {
"GHQ 0 (No distress)": "o",
"GHQ 1-3 (Mild)": "s",
"GHQ 4+ (Probable distress)": "^"
}
source_note = (
"Source: HSE 2022 (NHS England / UKDS, Study 9469). "
"Weighted using wt_int. Analytical sample N=4,735."
)
Appendix A: Data and Variable Construction¶
Appendix A.1: GHQ-12 Psychological Distress Score¶
In [12]:
# ================================================================
# Appendix A.1: GHQ-12 Psychological Distress Score
# ================================================================
# Validated 12-item psychological distress measure, administered
# via self-completion booklet (SC 13+) in HSE 2022.
# HSE provides pre-derived GHQ variables — retained as-is.
# Caseness threshold: score >= 4 (probable psychological distress).
#
# Derived variables used:
# ghq12scr — continuous score 0-12 (primary exposure)
# ghqg2 — grouped: 1=score 0, 2=score 1-3, 3=score 4+
# ghq — binary caseness: 0=score 0-3, 1=score 4+
# ================================================================
# Verify cleaning worked correctly
print("=== Value check after cleaning ===")
print(df['ghq12scr'].value_counts(dropna=False).sort_index())
# Range assertions
assert df['ghq12scr'].dropna().min() >= 0, "GHQ score below 0"
assert df['ghq12scr'].dropna().max() <= 12, "GHQ score above 12"
print("\nRange assertions passed: ghq12scr within [0, 12]")
# N valid / missing
print(f"\nN valid: {df['ghq12scr'].notna().sum():,}")
print(f"N missing: {df['ghq12scr'].isna().sum():,} ({df['ghq12scr'].isna().sum()/len(df)*100:.1f}%)")
# GHQ category labels for figures (used in C)
df['ghq_cat'] = df['ghqg2'].map({
1: '0 (No distress)',
2: '1-3 (Mild)',
3: '4+ (Probable distress)'
})
=== Value check after cleaning === ghq12scr 0.0 3488 1.0 769 2.0 424 3.0 313 4.0 223 5.0 161 6.0 147 7.0 116 8.0 107 9.0 87 10.0 90 11.0 88 12.0 99 NaN 1617 Name: count, dtype: int64 Range assertions passed: ghq12scr within [0, 12] N valid: 6,112 N missing: 1,617 (20.9%)
Appendix A.2: Physical Condition Count¶
In [14]:
# ================================================================
# Appendix A.2: Physical Condition Count
# ================================================================
# condlcnt is a pre-derived HSE variable counting the number of
# grouped longstanding condition categories reported by the
# respondent. Conditions are grouped into ICD-based categories.
# Range observed in cleaned adult sample: 0-6 condition categories.
#
# Derived variables used:
# condlcnt — count of condition categories (primary outcome)
# condlcnt2 — same, capped at 4+ (descriptives)
# limlast — limiting longstanding illness (robustness)
# multimorbid — binary: 1 if condlcnt >= 2 (robustness)
#
# Individual condition flags (descriptives only, not in regression):
# compexp1 — Diabetes
# compexp3 — IHD/stroke/angina
# compexp4 — Hypertension
# compexp6 — COPD/bronchitis
# compexp7 — Asthma
# compexp9 — Arthritis/rheumatism
# ================================================================
# Binary multimorbidity indicator
df['multimorbid'] = (df['condlcnt'] >= 2).astype(float)
df.loc[df['condlcnt'].isna(), 'multimorbid'] = np.nan
# Verify cleaning
print("=== Value check after cleaning ===")
print(df['condlcnt'].value_counts(dropna=False).sort_index())
# Range assertion
assert df['condlcnt'].dropna().min() >= 0, "condlcnt below 0"
print("\nRange assertion passed: condlcnt >= 0")
# N valid / missing
print(f"\nN valid: {df['condlcnt'].notna().sum():,}")
print(f"N missing: {df['condlcnt'].isna().sum():,} ({df['condlcnt'].isna().sum()/len(df)*100:.1f}%)")
# Mean and variance — overdispersion check for model choice justification
mean_cond = df['condlcnt'].mean()
var_cond = df['condlcnt'].var()
print(f"\n=== Overdispersion Check ===")
print(f"Mean: {mean_cond:.3f}")
print(f"Variance: {var_cond:.3f}")
print(f"Variance / Mean: {var_cond/mean_cond:.3f}")
print("(Ratio > 1 indicates overdispersion — supports negative binomial over Poisson)")
# Verify individual condition flag coding before reporting prevalence
print("\n=== Individual Condition Flag Value Counts ===")
flags = {
'compexp1': 'Diabetes',
'compexp3': 'IHD/stroke/angina',
'compexp4': 'Hypertension',
'compexp6': 'COPD/bronchitis',
'compexp7': 'Asthma',
'compexp9': 'Arthritis/rheumatism'
}
for var, label in flags.items():
print(f"\n{label} ({var})")
print(df[var].value_counts(dropna=False).sort_index())
# Prevalence conditional on valid coding confirmed above
print("\n=== Condition Prevalence (adults 16+) ===")
for var, label in flags.items():
n = (df[var] == 1).sum()
pct = n / df[var].notna().sum() * 100
print(f" {label:<25} n={n:,} ({pct:.1f}%)")
=== Value check after cleaning === condlcnt 0.0 4098 1.0 1839 2.0 1023 3.0 481 4.0 201 5.0 72 6.0 8 NaN 7 Name: count, dtype: int64 Range assertion passed: condlcnt >= 0 N valid: 7,722 N missing: 7 (0.1%) === Overdispersion Check === Mean: 0.847 Variance: 1.298 Variance / Mean: 1.533 (Ratio > 1 indicates overdispersion — supports negative binomial over Poisson) === Individual Condition Flag Value Counts === Diabetes (compexp1) compexp1 0.0 7324 1.0 398 NaN 7 Name: count, dtype: int64 IHD/stroke/angina (compexp3) compexp3 0.0 7553 1.0 169 NaN 7 Name: count, dtype: int64 Hypertension (compexp4) compexp4 0.0 7268 1.0 454 NaN 7 Name: count, dtype: int64 COPD/bronchitis (compexp6) compexp6 0.0 7546 1.0 176 NaN 7 Name: count, dtype: int64 Asthma (compexp7) compexp7 0.0 7351 1.0 371 NaN 7 Name: count, dtype: int64 Arthritis/rheumatism (compexp9) compexp9 0.0 7017 1.0 705 NaN 7 Name: count, dtype: int64 === Condition Prevalence (adults 16+) === Diabetes n=398 (5.2%) IHD/stroke/angina n=169 (2.2%) Hypertension n=454 (5.9%) COPD/bronchitis n=176 (2.3%) Asthma n=371 (4.8%) Arthritis/rheumatism n=705 (9.1%)
Appendix A.3: Area Deprivation — IMD Quintile¶
In [15]:
# ================================================================
# Appendix A.3: Area Deprivation — IMD Quintile (qimd19)
# ================================================================
# qimd19 is a pre-derived HSE variable based on the Index of
# Multiple Deprivation 2019 (IMD 2019), linked to respondents
# via their area of residence.
# IMD is an area-level measure of relative deprivation across
# seven domains: income, employment, education, health,
# crime, housing, and environment.
# Quintile 1 = least deprived, Quintile 5 = most deprived.
# This is an area-level, not individual-level, measure.
# Used as key moderator in regression and interaction term.
# ================================================================
# IMD quintile labels for figures
df['imd_label'] = df['qimd19'].map({
1: 'Q1 Least\ndeprived',
2: 'Q2',
3: 'Q3',
4: 'Q4',
5: 'Q5 Most\ndeprived'
})
# Verify cleaning
print("=== Value check after cleaning ===")
print(df['qimd19'].value_counts(dropna=False).sort_index())
# Range assertion
assert df['qimd19'].dropna().min() == 1, "IMD quintile below 1"
assert df['qimd19'].dropna().max() == 5, "IMD quintile above 5"
print("\nRange assertions passed: qimd19 within [1, 5]")
# N valid / missing
print(f"\nN valid: {df['qimd19'].notna().sum():,}")
print(f"N missing: {df['qimd19'].isna().sum():,} ({df['qimd19'].isna().sum()/len(df)*100:.1f}%)")
=== Value check after cleaning === qimd19 1.0 1796 2.0 1631 3.0 1535 4.0 1428 5.0 1339 Name: count, dtype: int64 Range assertions passed: qimd19 within [1, 5] N valid: 7,729 N missing: 0 (0.0%)
Appendix A.4: Demographics¶
In [19]:
# ================================================================
# Appendix A.4: Demographics
# ================================================================
# Three demographic variables used as controls in regression.
#
# Sex (Sex)
# 1 = Male (reference category)
# 2 = Female
# Source: household grid
#
# Age (ag16g10)
# 10-year bands for adults 16+
# 1=16-24, 2=25-34, 3=35-44, 4=45-54,
# 5=55-64, 6=65-74, 7=75+
# Reference category: 16-24
#
# Ethnicity (origin2)
# Grouped ethnic categories — labels verified from HSE
# data dictionary below
# Reference category: White (1)
# ================================================================
# Verify ethnicity labels from data dictionary
reader = pd.io.stata.StataReader('Data/stata/stata13/hse_2022_eul_v1.dta')
value_labels = reader.value_labels()
print("=== origin2 official value labels ===")
for k, v in sorted(value_labels['origin2'].items()):
print(f" {k}: {v}")
# Verify cleaning
print("\n=== Sex (Sex) ===")
print(df['Sex'].value_counts(dropna=False).sort_index())
print("\n=== Age bands (ag16g10) ===")
print(df['ag16g10'].value_counts(dropna=False).sort_index())
print("\n=== Ethnicity (origin2) ===")
print(df['origin2'].value_counts(dropna=False).sort_index())
# Missingness
print("\n=== Missingness ===")
for var in ['Sex', 'ag16g10', 'origin2']:
n_miss = df[var].isna().sum()
pct = n_miss / len(df) * 100
print(f" {var:<12} N missing: {n_miss:,} ({pct:.1f}%)")
=== origin2 official value labels === -9: Refused -8: Don't know -1: Not applicable 1: White 2: Black 3: Asian 4: Mixed/multiple ethnic background 5: Any other ethnic group === Sex (Sex) === Sex 1.0 3488 2.0 4241 Name: count, dtype: int64 === Age bands (ag16g10) === ag16g10 1.0 498 2.0 887 3.0 1149 4.0 1175 5.0 1380 6.0 1452 7.0 1188 Name: count, dtype: int64 === Ethnicity (origin2) === origin2 1.0 6750 2.0 209 3.0 498 4.0 134 5.0 70 NaN 68 Name: count, dtype: int64 === Missingness === Sex N missing: 0 (0.0%) ag16g10 N missing: 0 (0.0%) origin2 N missing: 68 (0.9%)
Appendix A.5 Socioeconomic Status¶
In [20]:
# ================================================================
# Appendix A.5: Socioeconomic Status
# ================================================================
# Three socioeconomic variables used as controls in regression,
# capturing different dimensions of socioeconomic position.
#
# Income (eqv5)
# Equivalised household income quintiles
# Accounts for household size and composition
# 1 = Lowest quintile (<= £19,180) (reference category)
# 2 = Second lowest (> £19,180 <= £27,705)
# 3 = Middle (> £27,705 <= £42,763)
# 4 = Second highest (> £42,763 <= £66,901)
# 5 = Highest (> £66,901)
# Note: -90 treated as missing (age of household member refused)
#
# Social class (nssec3)
# NS-SEC 3-category classification (individual)
# 1 = Managerial and professional (reference category)
# 2 = Intermediate occupations
# 3 = Routine and manual occupations
# Note: value 99 (Other) treated as missing
#
# Education (topqual3)
# Highest educational qualification
# 1 = NVQ4/NVQ5/Degree or equivalent (reference category)
# 2 = Higher education below degree
# 3 = NVQ3/GCE A Level equivalent
# 4 = NVQ2/GCE O Level equivalent
# 5 = NVQ1/CSE other grade equivalent
# 6 = Foreign/other
# 7 = No qualification
# ================================================================
# Verify labels from data dictionary
reader = pd.io.stata.StataReader('Data/stata/stata13/hse_2022_eul_v1.dta')
value_labels = reader.value_labels()
for var, label in [('eqv5','Income quintiles'),
('nssec3','NS-SEC'),
('topqual3','Education')]:
print(f"\n=== {label} ({var}) official labels ===")
if var in value_labels:
for k, v in sorted(value_labels[var].items()):
print(f" {k}: {v}")
# Value counts after cleaning
print("\n=== Income quintiles (eqv5) ===")
print(df['eqv5'].value_counts(dropna=False).sort_index())
print("\n=== NS-SEC (nssec3) ===")
print(df['nssec3'].value_counts(dropna=False).sort_index())
print("\n=== Education (topqual3) ===")
print(df['topqual3'].value_counts(dropna=False).sort_index())
# Missingness
print("\n=== Missingness ===")
for var in ['eqv5', 'nssec3', 'topqual3']:
n_miss = df[var].isna().sum()
pct = n_miss / len(df) * 100
print(f" {var:<12} N missing: {n_miss:,} ({pct:.1f}%)")
=== Income quintiles (eqv5) official labels === -90: Age of household member refused -9: Refused -8: Don't know -1: Not applicable 1: Lowest Quintile (<= £19,180) 2: Second lowest Quintile (> £19,180 <= £27,705) 3: Middle Quintile (> £27,705 <= £42,763) 4: Second highest Quintile (> £42,763 <= £66,901) 5: Highest Quintile (> £66,901) === NS-SEC (nssec3) official labels === -9: Refused -8: Don't know -1: Not applicable 1: Managerial and professional occupations 2: Intermediate occupations 3: Routine and manual occupations 99: Other === Education (topqual3) official labels === -9: Refused -8: Don't know -1: Not applicable 1: NVQ4/NVQ5/Degree or equiv 2: Higher ed below degree 3: NVQ3/GCE A Level equiv 4: NVQ2/GCE O Level equiv 5: NVQ1/CSE other grade equiv 6: Foreign/other 7: No qualification === Income quintiles (eqv5) === eqv5 1.0 1057 2.0 1273 3.0 1210 4.0 1278 5.0 1218 NaN 1693 Name: count, dtype: int64 === NS-SEC (nssec3) === nssec3 1.0 3079 2.0 1686 3.0 2442 NaN 522 Name: count, dtype: int64 === Education (topqual3) === topqual3 1.0 2628 2.0 882 3.0 1130 4.0 1357 5.0 234 6.0 68 7.0 1356 NaN 74 Name: count, dtype: int64 === Missingness === eqv5 N missing: 1,693 (21.9%) nssec3 N missing: 522 (6.8%) topqual3 N missing: 74 (1.0%)
Appendix A.6: Lifestyle¶
In [22]:
# ================================================================
# Appendix A.6: Lifestyle
# ================================================================
# Two lifestyle variables used as controls in regression,
# capturing health behaviours that are associated with both
# psychological distress and physical multimorbidity.
#
# Smoking (cigst1_19)
# Cigarette smoking status
# 1 = Never smoked (reference category)
# 2 = Used to smoke occasionally
# 3 = Used to smoke regularly
# 4 = Current smoker
#
# Alcohol (totalwug2_22)
# Alcohol consumption risk groups (2022 revised units)
# 0 = Non-drinker / not drunk in last 12 months
# (reference category)
# 1 = Lower risk (up to 14 units per week)
# 2 = Increased risk (over 14-35/50 units per week)
# 3 = Higher risk (more than 35/50 units per week)
# ================================================================
# Verify labels from data dictionary
reader = pd.io.stata.StataReader('Data/stata/stata13/hse_2022_eul_v1.dta')
value_labels = reader.value_labels()
for var, label in [('cigst1_19', 'Smoking'),
('totalwug2_22', 'Alcohol')]:
print(f"\n=== {label} ({var}) official labels ===")
if var in value_labels:
for k, v in sorted(value_labels[var].items()):
print(f" {k}: {v}")
# Value counts after cleaning
print("\n=== Smoking (cigst1_19) ===")
print(df['cigst1_19'].value_counts(dropna=False).sort_index())
print("\n=== Alcohol (totalwug2_22) ===")
print(df['totalwug2_22'].value_counts(dropna=False).sort_index())
# Missingness
print("\n=== Missingness ===")
for var in ['cigst1_19', 'totalwug2_22']:
n_miss = df[var].isna().sum()
pct = n_miss / len(df) * 100
print(f" {var:<15} N missing: {n_miss:,} ({pct:.1f}%)")
=== Smoking (cigst1_19) official labels === -9: Refused -8: Don't know -1: Not applicable 1: Never smoked cigarettes at all 2: Used to smoke cigarettes occasionally 3: Used to smoke cigarettes regularly 4: Current cigarette smoker === Alcohol (totalwug2_22) official labels === -9: Refused -8: Don't know -1: Not applicable 0: Non drinker/not in last 12 months 1: Lower risk (up to 14 units) 2: Increased risk (over 14-50/ over 14-35) 3: Higher risk (more than 50/35) === Smoking (cigst1_19) === cigst1_19 1.0 4130 2.0 570 3.0 2011 4.0 917 NaN 101 Name: count, dtype: int64 === Alcohol (totalwug2_22) === totalwug2_22 0.0 1371 1.0 4286 2.0 1438 3.0 362 NaN 272 Name: count, dtype: int64 === Missingness === cigst1_19 N missing: 101 (1.3%) totalwug2_22 N missing: 272 (3.5%)
Appendix A.7: Geography¶
In [23]:
# ================================================================
# Appendix A.7: Geography
# ================================================================
# One geographic variable used as a control in regression,
# capturing regional variation in health outcomes and
# healthcare provision across England.
#
# Region (GOR1)
# Government Office Region (9 regions)
# 1 = North East (reference category)
# 2 = North West
# 3 = Yorkshire and The Humber
# 4 = East Midlands
# 5 = West Midlands
# 6 = East of England
# 7 = London
# 8 = South East
# 9 = South West
# ================================================================
# Verify labels from data dictionary
reader = pd.io.stata.StataReader('Data/stata/stata13/hse_2022_eul_v1.dta')
value_labels = reader.value_labels()
print("=== Region (GOR1) official labels ===")
if 'GOR1' in value_labels:
for k, v in sorted(value_labels['GOR1'].items()):
print(f" {k}: {v}")
# Value counts after cleaning
print("\n=== Region (GOR1) ===")
print(df['GOR1'].value_counts(dropna=False).sort_index())
# Missingness
n_miss = df['GOR1'].isna().sum()
pct = n_miss / len(df) * 100
print(f"\nN valid: {df['GOR1'].notna().sum():,}")
print(f"N missing: {n_miss:,} ({pct:.1f}%)")
=== Region (GOR1) official labels === -9: Refused -8: Don't know -1: Not applicable 1: E12000001 North East 2: E12000002 North West 3: E12000003 Yorkshire and The Humber 4: E12000004 East Midlands 5: E12000005 West Midlands 6: E12000006 East of England 7: E12000007 London 8: E12000008 South East 9: E12000009 South West === Region (GOR1) === GOR1 1.0 804 2.0 911 3.0 780 4.0 569 5.0 720 6.0 849 7.0 826 8.0 1367 9.0 903 Name: count, dtype: int64 N valid: 7,729 N missing: 0 (0.0%)
Appendix A.8: BMI¶
In [24]:
# ================================================================
# Appendix A.8: BMI — Considered and Excluded
# ================================================================
# BMI was considered as a lifestyle control variable given its
# established association with both psychological distress and
# physical multimorbidity. However, it was excluded from the
# regression model due to prohibitive missingness.
#
# bmisrg3 — self-reported BMI grouped
# 1 = Underweight/normal
# 2 = Overweight
# 3 = Obese
# ================================================================
# Check missingness
print("=== BMI (bmisrg3) value counts ===")
print(df['bmisrg3'].value_counts(dropna=False).sort_index())
n_miss = df['bmisrg3'].isna().sum()
pct = n_miss / len(df) * 100
print(f"\nN valid: {df['bmisrg3'].notna().sum():,}")
print(f"N missing: {n_miss:,} ({pct:.1f}%)")
# Show impact on analytical sample if BMI were included
reg_vars_with_bmi = [
'ghq12scr', 'condlcnt', 'qimd19', 'Sex', 'ag16g10',
'origin2', 'eqv5', 'nssec3', 'topqual3',
'cigst1_19', 'totalwug2_22', 'GOR1', 'bmisrg3'
]
n_with_bmi = df.dropna(subset=reg_vars_with_bmi).shape[0]
reg_vars_without_bmi = [v for v in reg_vars_with_bmi if v != 'bmisrg3']
n_without_bmi = df.dropna(subset=reg_vars_without_bmi).shape[0]
print(f"\n=== Impact of Including BMI ===")
print(f"Analytical sample without BMI: {n_without_bmi:,}")
print(f"Analytical sample with BMI: {n_with_bmi:,}")
print(f"Sample loss from adding BMI: {n_without_bmi - n_with_bmi:,} ({(n_without_bmi - n_with_bmi)/n_without_bmi*100:.1f}%)")
print(f"\nBMI excluded on these grounds.")
=== BMI (bmisrg3) value counts === bmisrg3 1.0 1001 2.0 872 3.0 596 NaN 5260 Name: count, dtype: int64 N valid: 2,469 N missing: 5,260 (68.1%) === Impact of Including BMI === Analytical sample without BMI: 4,735 Analytical sample with BMI: 1,075 Sample loss from adding BMI: 3,660 (77.3%) BMI excluded on these grounds.
Appendix B: Sample Construction and Missingness¶
Appendix B.1: Sample Construction¶
In [26]:
# ================================================================
# Appendix B.1: Sample Construction
# ================================================================
reg_vars = [
'ghq12scr', 'condlcnt', 'qimd19', 'Sex', 'ag16g10',
'origin2', 'eqv5', 'nssec3', 'topqual3',
'cigst1_19', 'totalwug2_22', 'GOR1'
]
# Define analytical sample
df_analytic = df.dropna(subset=reg_vars).copy()
df_excluded = df[~df.index.isin(df_analytic.index)].copy()
# Sequential sample construction table
print("=== Table B.1: Sequential Sample Construction ===")
print(f"{'Step':<35} {'N Remaining':>12} {'N Lost':>10}")
print("-" * 60)
running = df.copy()
print(f"{'Full adult sample (16+)':<35} {len(running):>12} {'—':>10}")
var_order = [
('GHQ-12', 'ghq12scr'),
('Condition count', 'condlcnt'),
('IMD quintile', 'qimd19'),
('Income', 'eqv5'),
('NS-SEC', 'nssec3'),
('Education', 'topqual3'),
('Smoking', 'cigst1_19'),
('Alcohol', 'totalwug2_22'),
('Sex', 'Sex'),
('Age', 'ag16g10'),
('Ethnicity', 'origin2'),
('Region', 'GOR1'),
]
for label, var in var_order:
before = len(running)
running = running.dropna(subset=[var])
after = len(running)
lost = before - after
print(f"{'Drop missing: '+label:<35} {after:>12} {lost:>10}")
print("-" * 60)
print(f"{'Final analytical sample':<35} {len(df_analytic):>12}")
=== Table B.1: Sequential Sample Construction === Step N Remaining N Lost ------------------------------------------------------------ Full adult sample (16+) 7729 — Drop missing: GHQ-12 6112 1617 Drop missing: Condition count 6109 3 Drop missing: IMD quintile 6109 0 Drop missing: Income 5041 1068 Drop missing: NS-SEC 4841 200 Drop missing: Education 4830 11 Drop missing: Smoking 4817 13 Drop missing: Alcohol 4739 78 Drop missing: Sex 4739 0 Drop missing: Age 4739 0 Drop missing: Ethnicity 4735 4 Drop missing: Region 4735 0 ------------------------------------------------------------ Final analytical sample 4735
Appendix B.2: Missing Data Analysis¶
In [29]:
# ================================================================
# Appendix B.2: Missing Data Analysis
# ================================================================
compare_vars = [
('Age band', 'ag16g10'),
('Sex', 'Sex'),
('IMD quintile', 'qimd19'),
('Ethnicity', 'origin2'),
('Condition count', 'condlcnt'),
]
# --- Table B.2: Missingness Rates ---
print("=== Table B.2: Variable Missingness Rates (Adult Sample, N=7,729) ===")
print(f"{'Variable':<30} {'N Valid':>10} {'N Missing':>10} {'% Missing':>10}")
print("-" * 63)
miss_vars = [
('GHQ-12 (ghq12scr)', 'ghq12scr'),
('Condition count (condlcnt)', 'condlcnt'),
('IMD quintile (qimd19)', 'qimd19'),
('Income (eqv5)', 'eqv5'),
('NS-SEC (nssec3)', 'nssec3'),
('Education (topqual3)', 'topqual3'),
('Smoking (cigst1_19)', 'cigst1_19'),
('Alcohol (totalwug2_22)', 'totalwug2_22'),
('Sex', 'Sex'),
('Age (ag16g10)', 'ag16g10'),
('Ethnicity (origin2)', 'origin2'),
('Region (GOR1)', 'GOR1'),
]
for label, var in miss_vars:
n_valid = df[var].notna().sum()
n_miss = df[var].isna().sum()
pct = n_miss / len(df) * 100
print(f"{label:<30} {n_valid:>10,} {n_miss:>10,} {pct:>9.1f}%")
# --- Table B.3: Who is missing GHQ-12? ---
print("\n=== Table B.3: Characteristics by GHQ-12 Response Status ===")
print(f"{'Variable':<25} {'Has GHQ (n=6,112)':>20} {'Missing GHQ (n=1,617)':>22} {'Difference':>12}")
print("-" * 82)
ghq_present = df[df['ghq12scr'].notna()]
ghq_missing = df[df['ghq12scr'].isna()]
for label, var in compare_vars:
m_pres = ghq_present[var].mean()
m_miss = ghq_missing[var].mean()
diff = m_pres - m_miss
print(f"{label:<25} {m_pres:>20.3f} {m_miss:>22.3f} {diff:>+12.3f}")
# --- Table B.4: Who is missing income? ---
print("\n=== Table B.4: Characteristics by Income Response Status ===")
print(f"{'Variable':<25} {'Has Income (n=6,036)':>22} {'Missing Income (n=1,693)':>25} {'Difference':>12}")
print("-" * 87)
inc_present = df[df['eqv5'].notna()]
inc_missing = df[df['eqv5'].isna()]
for label, var in compare_vars:
m_pres = inc_present[var].mean()
m_miss = inc_missing[var].mean()
diff = m_pres - m_miss
print(f"{label:<25} {m_pres:>22.3f} {m_miss:>25.3f} {diff:>+12.3f}")
=== Table B.2: Variable Missingness Rates (Adult Sample, N=7,729) === Variable N Valid N Missing % Missing --------------------------------------------------------------- GHQ-12 (ghq12scr) 6,112 1,617 20.9% Condition count (condlcnt) 7,722 7 0.1% IMD quintile (qimd19) 7,729 0 0.0% Income (eqv5) 6,036 1,693 21.9% NS-SEC (nssec3) 7,207 522 6.8% Education (topqual3) 7,655 74 1.0% Smoking (cigst1_19) 7,628 101 1.3% Alcohol (totalwug2_22) 7,457 272 3.5% Sex 7,729 0 0.0% Age (ag16g10) 7,729 0 0.0% Ethnicity (origin2) 7,661 68 0.9% Region (GOR1) 7,729 0 0.0% === Table B.3: Characteristics by GHQ-12 Response Status === Variable Has GHQ (n=6,112) Missing GHQ (n=1,617) Difference ---------------------------------------------------------------------------------- Age band 4.502 4.224 +0.278 Sex 1.552 1.535 +0.017 IMD quintile 2.840 2.913 -0.072 Ethnicity 1.229 1.315 -0.087 Condition count 0.882 0.712 +0.170 === Table B.4: Characteristics by Income Response Status === Variable Has Income (n=6,036) Missing Income (n=1,693) Difference --------------------------------------------------------------------------------------- Age band 4.397 4.611 -0.214 Sex 1.545 1.561 -0.016 IMD quintile 2.814 3.002 -0.187 Ethnicity 1.225 1.323 -0.098 Condition count 0.842 0.863 -0.021
Appendix C: Descriptive Statistics and Figures¶
Appendix C.1: Sample Characteristics Table (Weighted)¶
In [76]:
# ================================================================
# Appendix C.1: Sample Characteristics Table
# ================================================================
def wmean(series, weights):
mask = series.notna() & weights.notna()
return np.average(series[mask], weights=weights[mask])
def wstd(series, weights):
mask = series.notna() & weights.notna()
avg = np.average(series[mask], weights=weights[mask])
variance = np.average((series[mask] - avg)**2, weights=weights[mask])
return np.sqrt(variance)
def wpct(series, value, weights):
mask = series.notna() & weights.notna()
w_val = weights[mask & (series == value)].sum()
w_tot = weights[mask].sum()
return w_val / w_tot * 100
w = df_analytic['wt_int']
print("=== Table C.1: Sample Characteristics (Analytical Sample, N=4,735) ===")
# --- Panel A: Continuous Variables ---
print("\nPanel A: Continuous Variables")
print(f"{'Variable':<30} {'W.Mean':>8} {'W.SD':>8} {'Min':>6} {'Max':>6}")
print("-" * 62)
for label, var in [('GHQ-12 score (0-12)', 'ghq12scr'),
('Condition count', 'condlcnt')]:
wm = wmean(df_analytic[var], w)
sd = wstd(df_analytic[var], w)
mn = df_analytic[var].min()
mx = df_analytic[var].max()
print(f"{label:<30} {wm:>8.3f} {sd:>8.3f} {mn:>6.0f} {mx:>6.0f}")
# --- Panel B: Categorical Variables ---
print("\nPanel B: Categorical Variables")
print(f"{'Variable':<35} {'Weighted %':>12}")
print("-" * 50)
# Sex
print("\nSex")
for val, label in [(1,'Male'), (2,'Female')]:
pct = wpct(df_analytic['Sex'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# Age group
print("\nAge group")
for val, label in [(1,'16-24'),(2,'25-34'),(3,'35-44'),(4,'45-54'),
(5,'55-64'),(6,'65-74'),(7,'75+')]:
pct = wpct(df_analytic['ag16g10'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# Ethnicity
print("\nEthnicity")
for val, label in [(1,'White'),(2,'Black'),(3,'Asian'),
(4,'Mixed/multiple'),(5,'Any other')]:
pct = wpct(df_analytic['origin2'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# IMD quintile
print("\nIMD quintile")
for val, label in [(1,'Q1 (Least deprived)'),(2,'Q2'),(3,'Q3'),
(4,'Q4'),(5,'Q5 (Most deprived)')]:
pct = wpct(df_analytic['qimd19'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# Income quintile
print("\nIncome quintile")
for val, label in [(1,'Q1 Lowest (<=£19,180)'),
(2,'Q2 (£19,181-£27,705)'),
(3,'Q3 (£27,706-£42,763)'),
(4,'Q4 (£42,764-£66,901)'),
(5,'Q5 Highest (>£66,901)')]:
pct = wpct(df_analytic['eqv5'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# NS-SEC
print("\nNS-SEC")
for val, label in [(1,'Managerial/professional'),
(2,'Intermediate'),
(3,'Routine/manual')]:
pct = wpct(df_analytic['nssec3'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# Education
print("\nEducation")
for val, label in [(1,'Degree or equivalent'),
(2,'Higher ed below degree'),
(3,'A Level equivalent'),
(4,'O Level equivalent'),
(5,'NVQ1/CSE'),
(6,'Foreign/other'),
(7,'No qualification')]:
pct = wpct(df_analytic['topqual3'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# Smoking
print("\nSmoking status")
for val, label in [(1,'Never smoked'),
(2,'Ex-occasional smoker'),
(3,'Ex-regular smoker'),
(4,'Current smoker')]:
pct = wpct(df_analytic['cigst1_19'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# Alcohol
print("\nAlcohol risk")
for val, label in [(0,'Non-drinker'),
(1,'Lower risk'),
(2,'Increased risk'),
(3,'Higher risk')]:
pct = wpct(df_analytic['totalwug2_22'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# GHQ caseness
print("\nGHQ-12 caseness")
for val, label in [(0,'Score 0-3 (no caseness)'),
(1,'Score 4+ (probable distress)')]:
pct = wpct(df_analytic['ghq'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
# Multimorbidity
print("\nMultimorbidity (2+ conditions)")
for val, label in [(0,'No (0-1 conditions)'),
(1,'Yes (2+ conditions)')]:
pct = wpct(df_analytic['multimorbid'], val, w)
print(f" {label:<33} {pct:>11.1f}%")
=== Table C.1: Sample Characteristics (Analytical Sample, N=4,735) === Panel A: Continuous Variables Variable W.Mean W.SD Min Max -------------------------------------------------------------- GHQ-12 score (0-12) 1.653 2.839 0 12 Condition count 0.759 1.093 0 6 Panel B: Categorical Variables Variable Weighted % -------------------------------------------------- Sex Male 48.6% Female 51.4% Age group 16-24 7.7% 25-34 17.7% 35-44 17.0% 45-54 16.4% 55-64 17.6% 65-74 13.6% 75+ 10.0% Ethnicity White 87.6% Black 3.1% Asian 6.2% Mixed/multiple 2.2% Any other 0.9% IMD quintile Q1 (Least deprived) 24.1% Q2 22.2% Q3 20.5% Q4 17.9% Q5 (Most deprived) 15.2% Income quintile Q1 Lowest (<=£19,180) 16.2% Q2 (£19,181-£27,705) 19.2% Q3 (£27,706-£42,763) 19.1% Q4 (£42,764-£66,901) 22.1% Q5 Highest (>£66,901) 23.5% NS-SEC Managerial/professional 46.4% Intermediate 21.6% Routine/manual 32.0% Education Degree or equivalent 41.2% Higher ed below degree 12.8% A Level equivalent 17.2% O Level equivalent 15.4% NVQ1/CSE 2.5% Foreign/other 0.5% No qualification 10.4% Smoking status Never smoked 53.8% Ex-occasional smoker 8.5% Ex-regular smoker 25.4% Current smoker 12.3% Alcohol risk Non-drinker 14.9% Lower risk 58.9% Increased risk 21.0% Higher risk 5.3% GHQ-12 caseness Score 0-3 (no caseness) 82.6% Score 4+ (probable distress) 17.4% Multimorbidity (2+ conditions) No (0-1 conditions) 79.7% Yes (2+ conditions) 20.3%
Appendix C.2: Weighted vs Unweighted Comparison¶
In [36]:
# ================================================================
# Appendix C.2: Weighted vs Unweighted Comparison
# ================================================================
print("=== Table C.2: Weighted vs Unweighted Comparison ===")
print(f"{'Variable':<35} {'Unweighted':>12} {'Weighted':>12} {'Difference':>12}")
print("-" * 74)
# Continuous variables
print("\nContinuous variables (mean)")
for label, var in [('GHQ-12 score', 'ghq12scr'),
('Condition count', 'condlcnt')]:
uw = df_analytic[var].mean()
wt = wmean(df_analytic[var], w)
print(f" {label:<33} {uw:>12.3f} {wt:>12.3f} {wt-uw:>+12.3f}")
# GHQ caseness
print("\nGHQ caseness (% score 4+)")
uw_ghq = (df_analytic['ghq'] == 1).sum() / df_analytic['ghq'].notna().sum() * 100
wt_ghq = wpct(df_analytic['ghq'], 1, w)
print(f" {'GHQ 4+ (probable distress)':<33} {uw_ghq:>11.1f}% {wt_ghq:>11.1f}% {wt_ghq-uw_ghq:>+11.1f}%")
# Multimorbidity
print("\nMultimorbidity (% 2+ conditions)")
uw_mm = (df_analytic['multimorbid'] == 1).sum() / df_analytic['multimorbid'].notna().sum() * 100
wt_mm = wpct(df_analytic['multimorbid'], 1, w)
print(f" {'2+ conditions':<33} {uw_mm:>11.1f}% {wt_mm:>11.1f}% {wt_mm-uw_mm:>+11.1f}%")
# IMD quintile distribution
print("\nIMD quintile (%)")
for val, label in [(1,'Q1 (Least deprived)'),(2,'Q2'),(3,'Q3'),
(4,'Q4'),(5,'Q5 (Most deprived)')]:
uw = (df_analytic['qimd19'] == val).sum() / len(df_analytic) * 100
wt = wpct(df_analytic['qimd19'], val, w)
print(f" {label:<33} {uw:>11.1f}% {wt:>11.1f}% {wt-uw:>+11.1f}%")
# Income quintile distribution
print("\nIncome quintile (%)")
for val, label in [(1,'Q1 Lowest'),(2,'Q2'),(3,'Q3'),
(4,'Q4'),(5,'Q5 Highest')]:
uw = (df_analytic['eqv5'] == val).sum() / len(df_analytic) * 100
wt = wpct(df_analytic['eqv5'], val, w)
print(f" {label:<33} {uw:>11.1f}% {wt:>11.1f}% {wt-uw:>+11.1f}%")
=== Table C.2: Weighted vs Unweighted Comparison === Variable Unweighted Weighted Difference -------------------------------------------------------------------------- Continuous variables (mean) GHQ-12 score 1.654 1.653 -0.000 Condition count 0.867 0.759 -0.108 GHQ caseness (% score 4+) GHQ 4+ (probable distress) 17.2% 17.4% +0.1% Multimorbidity (% 2+ conditions) 2+ conditions 23.8% 20.3% -3.5% IMD quintile (%) Q1 (Least deprived) 24.7% 24.1% -0.6% Q2 22.2% 22.2% +0.1% Q3 19.4% 20.5% +1.1% Q4 17.7% 17.9% +0.2% Q5 (Most deprived) 16.0% 15.2% -0.8% Income quintile (%) Q1 Lowest 15.8% 16.2% +0.4% Q2 21.5% 19.2% -2.3% Q3 20.1% 19.1% -1.0% Q4 21.8% 22.1% +0.3% Q5 Highest 20.9% 23.5% +2.6%
Appendix C.3: GHQ-12 Distribution¶
In [51]:
# ================================================================
# Appendix C.3: GHQ-12 Distribution
# ================================================================
from scipy.stats import skew
ghq_vals = df_analytic['ghq12scr'].dropna()
w_ghq = df_analytic.loc[ghq_vals.index, 'wt_int']
# Summary statistics
print("=== GHQ-12 Summary Statistics (Analytical Sample) ===")
print(f"N valid: {len(ghq_vals):,}")
print(f"Weighted mean: {wmean(df_analytic['ghq12scr'], w):.3f}")
print(f"Median: {ghq_vals.median():.0f}")
print(f"SD (unweighted): {ghq_vals.std():.3f}")
print(f"Skewness: {skew(ghq_vals):.3f}")
print(f"Proportion at 0: {(ghq_vals == 0).mean() * 100:.1f}%")
print(f"Proportion at 12: {(ghq_vals == 12).mean() * 100:.1f}%")
print(f"Weighted caseness: {wpct(df_analytic['ghq'], 1, w):.1f}%")
# Weighted histogram
fig, ax = plt.subplots(figsize=(9, 5))
ax.hist(ghq_vals, bins=13, range=(-0.5, 12.5),
weights=w_ghq,
color='#4393c3', edgecolor='white', linewidth=0.5, alpha=0.85)
ax.axvline(x=3.5, color='#d6604d', linestyle='--', linewidth=1.5,
label='Caseness threshold (score ≥ 4)')
ax.set_xlabel('GHQ-12 Score (0 = no distress, 12 = severe distress)', fontsize=11)
ax.set_ylabel('Weighted Frequency', fontsize=11)
ax.set_title('Distribution of GHQ-12 Psychological Distress Scores\nAdults Aged 16+ in England (HSE 2022)',
fontsize=11, fontweight='bold')
ax.set_xticks(range(0, 13))
ax.legend(fontsize=10)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.grid(axis='y', alpha=0.3, linestyle='--')
fig.text(0.12, 0.01,
'Source: Health Survey for England 2022 (NHS England / UK Data Service, Study 9469). Weighted using wt_int.',
fontsize=8.5, color='#555555')
plt.tight_layout()
plt.savefig('figure_c3_ghq_distribution.png', dpi=200, bbox_inches='tight')
plt.show()
print("\nFigure C.3 saved.")
=== GHQ-12 Summary Statistics (Analytical Sample) === N valid: 4,735 Weighted mean: 1.653 Median: 0 SD (unweighted): 2.867 Skewness: 2.038 Proportion at 0: 58.1% Proportion at 12: 1.4% Weighted caseness: 17.4%
Figure C.3 saved.
Appendix C.4: Condition Count Distribution¶
In [52]:
# ================================================================
# Appendix C.4: Condition Count Distribution
# ================================================================
from scipy.stats import skew
cond_vals = df_analytic['condlcnt'].dropna()
w_cond = df_analytic.loc[cond_vals.index, 'wt_int']
# Summary statistics
print("=== Condition Count Summary Statistics (Analytical Sample) ===")
print(f"N valid: {len(cond_vals):,}")
print(f"Median: {cond_vals.median():.0f}")
print(f"SD (unweighted): {cond_vals.std():.3f}")
print(f"Skewness: {skew(cond_vals):.3f}")
print(f"Proportion at 0: {(cond_vals == 0).mean() * 100:.1f}%")
print(f"Proportion at 6: {(cond_vals == 6).mean() * 100:.1f}%")
print(f"Weighted mean: {wmean(df_analytic['condlcnt'], w):.3f}")
print(f"Variance: {cond_vals.var():.3f}")
print(f"Variance/Mean: {cond_vals.var() / cond_vals.mean():.3f}")
print(f"Weighted multimorbidity (2+): {wpct(df_analytic['multimorbid'], 1, w):.1f}%")
# Weighted histogram
fig, ax = plt.subplots(figsize=(9, 5))
ax.hist(cond_vals, bins=7, range=(-0.5, 6.5),
weights=w_cond,
color='#4393c3', edgecolor='white', linewidth=0.5, alpha=0.85)
ax.axvline(x=1.5, color='#d6604d', linestyle='--', linewidth=1.5,
label='Multimorbidity threshold (2+ conditions)')
ax.set_xlabel('Number of Physical Longstanding Condition Categories', fontsize=11)
ax.set_ylabel('Weighted Frequency', fontsize=11)
ax.set_title('Distribution of Physical Condition Count\nAdults Aged 16+ in England (HSE 2022)',
fontsize=11, fontweight='bold')
ax.set_xticks(range(0, 7))
ax.legend(fontsize=10)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.grid(axis='y', alpha=0.3, linestyle='--')
fig.text(0.12, 0.01,
'Source: Health Survey for England 2022 (NHS England / UK Data Service, Study 9469). Weighted using wt_int.',
fontsize=8.5, color='#555555')
plt.tight_layout()
plt.savefig('figure_c4_condition_distribution.png', dpi=200, bbox_inches='tight')
plt.show()
print("\nFigure C.4 saved.")
=== Condition Count Summary Statistics (Analytical Sample) === N valid: 4,735 Median: 0 SD (unweighted): 1.155 Skewness: 1.416 Proportion at 0: 52.5% Proportion at 6: 0.1% Weighted mean: 0.759 Variance: 1.333 Variance/Mean: 1.538 Weighted multimorbidity (2+): 20.3%
Figure C.4 saved.
Appendix C.5: Condition Prevalence by GHQ Group¶
In [ ]:
# ================================================================
# Appendix C.5: Condition Prevalence by GHQ Group
# ================================================================
df_analytic['ghq_cat'] = df_analytic['ghqg2'].map({
1: '0 (No distress)',
2: '1-3 (Mild)',
3: '4+ (Probable distress)'
})
flags = {
'compexp1': 'Diabetes',
'compexp3': 'IHD/stroke/angina',
'compexp4': 'Hypertension',
'compexp6': 'COPD/bronchitis',
'compexp7': 'Asthma',
'compexp9': 'Arthritis/rheumatism'
}
print("=== Table C.5: Condition Prevalence by GHQ-12 Group (Weighted %) ===")
print(f"{'Condition':<25} {'GHQ 0':>10} {'GHQ 1-3':>10} {'GHQ 4+':>10} {'Total':>10}")
print("-" * 60)
for var, label in flags.items():
row = []
for ghq_val in [1, 2, 3]:
subset = df_analytic[df_analytic['ghqg2'] == ghq_val]
w_sub = subset['wt_int']
pct = wpct(subset[var], 1, w_sub)
row.append(pct)
total = wpct(df_analytic[var], 1, w)
print(f"{label:<25} {row[0]:>9.1f}% {row[1]:>9.1f}% {row[2]:>9.1f}% {total:>9.1f}%")
# Weighted mean condition count by GHQ group
print("\n=== Mean Condition Count by GHQ Group (Weighted) ===")
print(f"{'GHQ Group':<25} {'N':>8} {'W.Mean condlcnt':>17} {'W.% Multimorbid':>17}")
print("-" * 70)
for ghq_val, label in [(1,'0 (No distress)'),
(2,'1-3 (Mild)'),
(3,'4+ (Probable distress)')]:
subset = df_analytic[df_analytic['ghqg2'] == ghq_val]
w_sub = subset['wt_int']
n = len(subset)
wm = wmean(subset['condlcnt'], w_sub)
mm = wpct(subset['multimorbid'], 1, w_sub)
print(f"{label:<25} {n:>8,} {wm:>17.3f} {mm:>16.1f}%")
=== Table C.5: Condition Prevalence by GHQ-12 Group (Weighted %) === Condition GHQ 0 GHQ 1-3 GHQ 4+ Total ------------------------------------------------------------ Diabetes 4.3% 3.3% 5.7% 4.3% IHD/stroke/angina 1.5% 1.1% 2.2% 1.5% Hypertension 5.9% 3.8% 4.7% 5.2% COPD/bronchitis 1.3% 1.8% 2.3% 1.6% Asthma 3.8% 6.1% 7.2% 5.0% Arthritis/rheumatism 6.2% 6.8% 12.5% 7.5% === Mean Condition Count by GHQ Group (Weighted) === GHQ Group N W.Mean condlcnt W.% Multimorbid ---------------------------------------------------------------------- 0 (No distress) 2,753 0.594 15.3% 1-3 (Mild) 1,166 0.783 21.0% 4+ (Probable distress) 816 1.270 36.1%
In [80]:
print("=== Table C.5c: Mean Condition Count and Multimorbidity by IMD Quintile (Weighted) ===")
print(f"{'IMD Quintile':<25} {'N':>8} {'W.Mean condlcnt':>17} {'W.% Multimorbid':>17}")
print("-" * 70)
for q, lab in [(1,'Q1 (Least deprived)'),(2,'Q2'),(3,'Q3'),
(4,'Q4'),(5,'Q5 (Most deprived)')]:
subset = df_analytic[df_analytic['qimd19'] == q]
w_sub = subset['wt_int']
n = len(subset)
wm = wmean(subset['condlcnt'], w_sub)
mm = wpct(subset['multimorbid'], 1, w_sub)
print(f"{lab:<25} {n:>8,} {wm:>17.3f} {mm:>16.1f}%")
=== Table C.5c: Mean Condition Count and Multimorbidity by IMD Quintile (Weighted) === IMD Quintile N W.Mean condlcnt W.% Multimorbid ---------------------------------------------------------------------- Q1 (Least deprived) 1,168 0.623 15.6% Q2 1,051 0.759 21.1% Q3 920 0.758 20.1% Q4 837 0.744 19.8% Q5 (Most deprived) 759 0.991 27.6%
Appendix C.6: Figures¶
In [82]:
# ================================================================
# Appendix C.6: Descriptive Figures
# ================================================================
import textwrap
source_note_fig1 = (
"Source: NHS England (2022) Health Survey for England, 2022 [data collection]. "
"UK Data Service, SN 9469. Weighted estimates using HSE 2022 interview weights (wt_int). "
"GHQ-12 caseness threshold: score ≥ 4 indicates probable psychological distress "
"(Goldberg et al., 1997). IMD: Index of Multiple Deprivation 2019 (MHCLG, 2019)."
)
source_note_c61 = (
"Source: NHS England (2022) Health Survey for England, 2022 [data collection]. "
"UK Data Service, SN 9469. Weighted using wt_int. "
"Conditions derived from HSE compexp variables (pre-coded ICD-based condition flags)."
)
source_note_c62 = (
"Source: NHS England (2022) Health Survey for England, 2022 [data collection]. "
"UK Data Service, SN 9469. Weighted using wt_int. "
"IMD: Index of Multiple Deprivation 2019 (MHCLG, 2019)."
)
GHQ_COLORS = {
"GHQ 0 (No distress)": "#2C7FB8",
"GHQ 1-3 (Mild)": "#7FCDBB",
"GHQ 4+ (Probable distress)": "#D95F0E"
}
GHQ_MARKERS = {
"GHQ 0 (No distress)": "o",
"GHQ 1-3 (Mild)": "s",
"GHQ 4+ (Probable distress)": "^"
}
# ================================================================
# Figure 1 (Main Report Figure)
# ================================================================
line_data = []
for ghq_val, ghq_lab in [
(1, "GHQ 0 (No distress)"),
(2, "GHQ 1-3 (Mild)"),
(3, "GHQ 4+ (Probable distress)")
]:
for imd_val in range(1, 6):
subset = df_analytic[
(df_analytic["ghqg2"] == ghq_val) &
(df_analytic["qimd19"] == imd_val)
]
line_data.append({
"ghq": ghq_lab,
"imd": imd_val,
"mean_cond": wmean(subset["condlcnt"], subset["wt_int"]),
"n": len(subset)
})
line_df = pd.DataFrame(line_data)
fig, ax = plt.subplots(figsize=(9.5, 5.8))
for ghq_lab in GHQ_COLORS.keys():
sub = line_df[line_df["ghq"] == ghq_lab]
ax.plot(
sub["imd"], sub["mean_cond"],
color=GHQ_COLORS[ghq_lab],
marker=GHQ_MARKERS[ghq_lab],
linewidth=2.6, markersize=7.5,
label=ghq_lab
)
ghq4 = line_df[line_df["ghq"] == "GHQ 4+ (Probable distress)"]
for _, row in ghq4.iterrows():
ax.text(
row["imd"], row["mean_cond"] + 0.045,
f"{row['mean_cond']:.2f}",
ha="center", va="bottom", fontsize=9,
color=GHQ_COLORS["GHQ 4+ (Probable distress)"]
)
ax.set_xticks(range(1, 6))
ax.set_xticklabels(["Q1\nLeast deprived", "Q2", "Q3", "Q4", "Q5\nMost deprived"])
ax.set_ylim(0.45, 1.85)
ax.set_xlabel("Area Deprivation Quintile (IMD 2019)")
ax.set_ylabel("Weighted Mean Physical Condition Count")
ax.set_title(
"Figure 1: Mean Physical Condition Count by Psychological Distress and Area Deprivation\n"
"Adults Aged 16+ in England, HSE 2022",
fontweight="bold", pad=12
)
ax.legend(title="GHQ-12 Distress Level", loc="upper left",
frameon=True, framealpha=0.95)
ax.grid(axis="y", linestyle="--", alpha=0.25)
ax.grid(axis="x", linestyle="--", alpha=0.10)
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
wrapped_fig1 = "\n".join(textwrap.wrap(source_note_fig1, width=110))
fig.text(0.12, 0.01, wrapped_fig1, fontsize=8, color="#555555")
plt.tight_layout(rect=[0, 0.10, 1, 1])
plt.savefig("figure1.png", dpi=300, bbox_inches="tight")
plt.show()
print("Figure 1 saved.")
# ================================================================
# Figure C.6.1 (Appendix)
# Physical Condition Prevalence Heatmap by GHQ Group
# ================================================================
flags = {
"compexp1": "Diabetes",
"compexp3": "IHD/stroke/angina",
"compexp4": "Hypertension",
"compexp6": "COPD/bronchitis",
"compexp7": "Asthma",
"compexp9": "Arthritis/rheumatism"
}
conditions = list(flags.keys())
cond_labels = list(flags.values())
ghq_groups = [1, 2, 3]
ghq_labels = [
"GHQ 0\nNo distress",
"GHQ 1-3\nMild",
"GHQ 4+\nProbable distress"
]
heatmap_data = np.zeros((len(conditions), len(ghq_groups)))
for i, var in enumerate(conditions):
for j, ghq_val in enumerate(ghq_groups):
subset = df_analytic[df_analytic["ghqg2"] == ghq_val]
heatmap_data[i, j] = wpct(subset[var], 1, subset["wt_int"])
fig, ax = plt.subplots(figsize=(8.2, 5.8))
im = ax.imshow(heatmap_data, cmap="YlOrBr", aspect="auto",
vmin=0, vmax=np.nanmax(heatmap_data))
ax.set_xticks(np.arange(len(ghq_labels)))
ax.set_xticklabels(ghq_labels)
ax.set_yticks(np.arange(len(cond_labels)))
ax.set_yticklabels(cond_labels)
for i in range(len(conditions)):
for j in range(len(ghq_groups)):
value = heatmap_data[i, j]
ax.text(j, i, f"{value:.1f}%",
ha="center", va="center", fontsize=10,
color="white" if value >= 10 else "black")
cbar = fig.colorbar(im, ax=ax, fraction=0.045, pad=0.04)
cbar.set_label("Weighted Prevalence (%)")
ax.set_title(
"Figure C.6.1: Physical Condition Prevalence by Psychological Distress Level\n"
"Adults Aged 16+ in England, HSE 2022",
fontweight="bold", pad=12
)
wrapped_c61 = "\n".join(textwrap.wrap(source_note_c61, width=110))
fig.text(0.12, 0.01, wrapped_c61, fontsize=8, color="#555555")
plt.tight_layout(rect=[0, 0.10, 1, 1])
plt.savefig("figure_c6_1.png", dpi=300, bbox_inches="tight")
plt.show()
print("Figure C.6.1 saved.")
# ================================================================
# Figure C.6.2 (Appendix)
# Prevalence of Probable Psychological Distress by IMD Quintile
# ================================================================
imd_ghq = []
for imd_val in range(1, 6):
subset = df_analytic[df_analytic["qimd19"] == imd_val]
imd_ghq.append({
"imd": imd_val,
"pct": wpct(subset["ghq"], 1, subset["wt_int"])
})
imd_ghq_df = pd.DataFrame(imd_ghq)
fig, ax = plt.subplots(figsize=(8.5, 5))
bars = ax.bar(
imd_ghq_df["imd"], imd_ghq_df["pct"],
color="#2C7FB8", edgecolor="white",
linewidth=0.6, alpha=0.90
)
for bar, pct in zip(bars, imd_ghq_df["pct"]):
ax.text(
bar.get_x() + bar.get_width() / 2,
bar.get_height() + 0.4,
f"{pct:.1f}%",
ha="center", va="bottom", fontsize=10
)
ax.set_xticks(range(1, 6))
ax.set_xticklabels(["Q1\nLeast deprived", "Q2", "Q3", "Q4", "Q5\nMost deprived"])
ax.set_ylim(0, 26)
ax.set_xlabel("Area Deprivation Quintile (IMD 2019)")
ax.set_ylabel("Weighted % with GHQ-12 Score ≥ 4")
ax.set_title(
"Figure C.6.2: Prevalence of Probable Psychological Distress by Area Deprivation\n"
"Adults Aged 16+ in England, HSE 2022",
fontweight="bold", pad=12
)
ax.grid(axis="y", linestyle="--", alpha=0.25)
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
wrapped_c62 = "\n".join(textwrap.wrap(source_note_c62, width=110))
fig.text(0.12, 0.01, wrapped_c62, fontsize=8, color="#555555")
plt.tight_layout(rect=[0, 0.10, 1, 1])
plt.savefig("figure_c6_2.png", dpi=300, bbox_inches="tight")
plt.show()
print("Figure C.6.2 saved.")
Figure 1 saved.
Figure C.6.1 saved.
Figure C.6.2 saved.
In [78]:
# ================================================================
# Figure C.6.3: Geographic Analysis
# Mean Condition Count and Distress Prevalence by Region
# ================================================================
gor_labels = {
1: 'North East',
2: 'North West',
3: 'Yorkshire',
4: 'East Midlands',
5: 'West Midlands',
6: 'East of England',
7: 'London',
8: 'South East',
9: 'South West'
}
regional_data = []
for gor_val, gor_lab in gor_labels.items():
subset = df_analytic[df_analytic['GOR1'] == gor_val]
w_sub = subset['wt_int']
regional_data.append({
'region': gor_lab,
'mean_cond': wmean(subset['condlcnt'], w_sub),
'pct_ghq': wpct(subset['ghq'], 1, w_sub),
'n': len(subset)
})
reg_df = pd.DataFrame(regional_data).sort_values('mean_cond', ascending=True)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5.5))
# Panel A: Mean condition count
bars1 = ax1.barh(reg_df['region'], reg_df['mean_cond'],
color='#2C7FB8', edgecolor='white',
linewidth=0.5, alpha=0.88)
for bar, val in zip(bars1, reg_df['mean_cond']):
ax1.text(bar.get_width() + 0.01, bar.get_y() + bar.get_height()/2,
f'{val:.2f}', va='center', fontsize=9)
ax1.set_xlabel('Weighted Mean Physical Condition Count', fontsize=10)
ax1.set_title('A. Mean Physical Condition Count\nby Region', fontsize=10,
fontweight='bold')
ax1.set_xlim(0, 1.1)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
ax1.grid(axis='x', linestyle='--', alpha=0.25)
# Panel B: GHQ caseness — sorted same order as Panel A
reg_df2 = reg_df.copy()
bars2 = ax2.barh(reg_df2['region'], reg_df2['pct_ghq'],
color='#D95F0E', edgecolor='white',
linewidth=0.5, alpha=0.88)
for bar, val in zip(bars2, reg_df2['pct_ghq']):
ax2.text(bar.get_width() + 0.3, bar.get_y() + bar.get_height()/2,
f'{val:.1f}%', va='center', fontsize=9)
ax2.set_xlabel('Weighted % with GHQ-12 Score ≥ 4', fontsize=10)
ax2.set_title('B. Probable Psychological Distress\nby Region', fontsize=10,
fontweight='bold')
ax2.set_xlim(0, 28)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)
ax2.grid(axis='x', linestyle='--', alpha=0.25)
source_note_c63 = (
"Source: NHS England (2022) Health Survey for England, 2022 [data collection]. "
"UK Data Service, SN 9469. Weighted using wt_int. "
"Regions sorted by ascending mean physical condition count."
)
wrapped_c63 = "\n".join(textwrap.wrap(source_note_c63, width=120))
fig.text(0.12, 0.01, wrapped_c63, fontsize=7.5, color="#555555")
fig.suptitle(
'Figure C.6.3: Geographic Variation in Physical Multimorbidity and Psychological Distress\n'
'Adults Aged 16+ in England, HSE 2022',
fontsize=11, fontweight='bold', y=1.02
)
plt.tight_layout()
plt.savefig('figure_c6_3.png', dpi=300, bbox_inches='tight')
plt.show()
print("Figure C.6.3 saved.")
Figure C.6.3 saved.
Appendix D: Regression Results¶
Appendix D.1: Regression Table¶
In [66]:
# ================================================================
# Appendix D.1: Full Regression Results
# ================================================================
def fit_nb(xvars):
X = pd.DataFrame(
sm.add_constant(r[xvars].astype(float).values, has_constant='add'),
columns=['const'] + xvars
)
result = NegativeBinomial(r['condlcnt'].values, X.values).fit(
method='bfgs', maxiter=500, disp=False
)
param_names = ['const'] + xvars + ['alpha']
result.params = pd.Series(result.params, index=param_names)
result.bse = pd.Series(result.bse, index=param_names)
result.pvalues = pd.Series(result.pvalues, index=param_names)
return result
def get_bic(mod):
return -2 * mod.llf + mod.params.shape[0] * np.log(mod.nobs)
print("Fitting models...")
m1 = fit_nb(['ghq12scr'])
m2 = fit_nb(['ghq12scr'] + demo)
m3 = fit_nb(['ghq12scr'] + demo + ses + life + geo + imd)
m4 = fit_nb(['ghq12scr'] + demo + ses + life + geo + imd + inter)
print("Done.")
# --- Table D.1a: Model Fit Statistics ---
print("\n=== Table D.1a: Model Fit Statistics ===")
print(f"{'Model':<45} {'N':>8} {'Log-lik':>12} {'AIC':>10} {'BIC':>10} {'Alpha':>10}")
print("-" * 98)
for name, mod in [
('Model 1: GHQ only', m1),
('Model 2: + Demographics', m2),
('Model 3: + Full controls', m3),
('Model 4: + GHQ x IMD interaction', m4),
]:
bic = get_bic(mod)
print(f"{name:<45} {int(mod.nobs):>8} "
f"{mod.llf:>12.3f} {mod.aic:>10.3f} "
f"{bic:>10.3f} {mod.params['alpha']:>10.4f}")
# --- Table D.1b: GHQ-12 Coefficient Progression ---
print("\n=== Table D.1b: GHQ-12 Coefficient Across Progressive Models ===")
print("Note: Unweighted regression. Weighted descriptives in Appendix C.")
print(" In Model 4, ghq12scr = GHQ slope in IMD Q1 only (reference group).")
print(" Model 3 = main average association model.")
print(f"\n{'Model':<45} {'Coef':>8} {'IRR':>8} {'SE':>8} {'95% CI':>16} {'p':>8} {'':>5}")
print("-" * 100)
for name, mod in [
('Model 1: GHQ only', m1),
('Model 2: + Demographics', m2),
('Model 3: + Full controls', m3),
('Model 4: + GHQ x IMD (ref: Q1)', m4),
]:
c = mod.params['ghq12scr']
se = mod.bse['ghq12scr']
p = mod.pvalues['ghq12scr']
irr = np.exp(c)
ci_lo = np.exp(c - 1.96*se)
ci_hi = np.exp(c + 1.96*se)
stars = '***' if p<0.001 else '**' if p<0.01 else '*' if p<0.05 else ''
print(f"{name:<45} {c:>8.4f} {irr:>8.4f} {se:>8.4f} "
f"[{ci_lo:.3f}, {ci_hi:.3f}] {p:>8.4f} {stars:>5}")
# --- Full coefficient tables for all four models ---
var_labels = (
[('GHQ-12 score', 'ghq12scr')] +
[('Female (ref: Male)', 'female')] +
[(f'Age {l} (ref: 16-24)', f'age_{i}')
for i, l in zip(range(2,8),
['25-34','35-44','45-54','55-64','65-74','75+'])] +
[(f'Ethnicity: {l} (ref: White)', f'eth_{i}')
for i, l in zip(range(2,6),
['Black','Asian','Mixed','Other'])] +
[(f'Income {l} (ref: Q1 lowest)', f'inc_{i}')
for i, l in zip(range(2,6), ['Q2','Q3','Q4','Q5'])] +
[(f'NS-SEC: {l} (ref: Managerial)', f'nssec_{i}')
for i, l in zip(range(2,4),
['Intermediate','Routine/manual'])] +
[(f'Education: {l} (ref: Degree)', f'edu_{i}')
for i, l in zip(range(2,8),
['HE below degree','A Level','O Level',
'NVQ1/CSE','Foreign/other','No qual'])] +
[(f'Smoking: {l} (ref: Never)', f'smk_{i}')
for i, l in zip(range(2,5),
['Ex-occasional','Ex-regular','Current'])] +
[(f'Alcohol: {l} (ref: Lower risk)', f'alc_{i}')
for i, l in zip([0,2,3],
['Non-drinker','Increased risk','Higher risk'])] +
[(f'Region: {l} (ref: London)', f'gor_{i}')
for i, l in zip([1,2,3,4,5,6,8,9],
['North East','North West','Yorkshire',
'East Midlands','West Midlands',
'East of England','South East','South West'])] +
[(f'IMD {l} (ref: Q1 least deprived)', f'imd_{i}')
for i, l in zip(range(2,6), ['Q2','Q3','Q4','Q5'])]
)
for mod_name, mod, extra_vars in [
('Model 1: GHQ only', m1, []),
('Model 2: + Demographics', m2, []),
('Model 3: + Full controls', m3, []),
('Model 4: + GHQ x IMD interaction', m4,
[(f'GHQ x IMD {l} (ref: Q1)', f'ghq_x_imd{i}')
for i, l in zip(range(2,6), ['Q2','Q3','Q4','Q5'])]),
]:
print(f"\n=== Full Coefficients: {mod_name} ===")
print(f"{'Variable':<45} {'Coef':>8} {'IRR':>8} {'SE':>8} {'95% CI':>16} {'p':>8} {'':>5}")
print("-" * 100)
all_labels = var_labels + extra_vars
for label, var in all_labels:
if var not in mod.params.index:
continue
c = mod.params[var]
se = mod.bse[var]
p = mod.pvalues[var]
irr = np.exp(c)
ci_lo = np.exp(c - 1.96*se)
ci_hi = np.exp(c + 1.96*se)
stars = '***' if p<0.001 else '**' if p<0.01 else '*' if p<0.05 else ''
print(f"{label:<45} {c:>8.4f} {irr:>8.4f} {se:>8.4f} "
f"[{ci_lo:.3f}, {ci_hi:.3f}] {p:>8.4f} {stars:>5}")
bic = get_bic(mod)
print(f"\n N={int(mod.nobs):,} Log-lik={mod.llf:.3f} "
f"AIC={mod.aic:.3f} BIC={bic:.3f} "
f"Alpha={mod.params['alpha']:.4f}")
Fitting models...
Done.
=== Table D.1a: Model Fit Statistics ===
Model N Log-lik AIC BIC Alpha
--------------------------------------------------------------------------------------------------
Model 1: GHQ only 4735 -5976.887 11959.773 11979.161 0.6017
Model 2: + Demographics 4735 -5748.746 11525.492 11615.971 0.3884
Model 3: + Full controls 4735 -5658.049 11404.098 11688.459 0.3041
Model 4: + GHQ x IMD interaction 4735 -5657.445 11410.889 11721.101 0.3035
=== Table D.1b: GHQ-12 Coefficient Across Progressive Models ===
Note: Unweighted regression. Weighted descriptives in Appendix C.
In Model 4, ghq12scr = GHQ slope in IMD Q1 only (reference group).
Model 3 = main average association model.
Model Coef IRR SE 95% CI p
----------------------------------------------------------------------------------------------------
Model 1: GHQ only 0.0923 1.0967 0.0059 [1.084, 1.109] 0.0000 ***
Model 2: + Demographics 0.1043 1.1100 0.0056 [1.098, 1.122] 0.0000 ***
Model 3: + Full controls 0.0921 1.0965 0.0055 [1.085, 1.108] 0.0000 ***
Model 4: + GHQ x IMD (ref: Q1) 0.0936 1.0981 0.0132 [1.070, 1.127] 0.0000 ***
=== Full Coefficients: Model 1: GHQ only ===
Variable Coef IRR SE 95% CI p
----------------------------------------------------------------------------------------------------
GHQ-12 score 0.0923 1.0967 0.0059 [1.084, 1.109] 0.0000 ***
N=4,735 Log-lik=-5976.887 AIC=11959.773 BIC=11979.161 Alpha=0.6017
=== Full Coefficients: Model 2: + Demographics ===
Variable Coef IRR SE 95% CI p
----------------------------------------------------------------------------------------------------
GHQ-12 score 0.1043 1.1100 0.0056 [1.098, 1.122] 0.0000 ***
Female (ref: Male) 0.0442 1.0451 0.0374 [0.971, 1.125] 0.2380
Age 25-34 (ref: 16-24) -0.0111 0.9890 0.1394 [0.753, 1.300] 0.9368
Age 35-44 (ref: 16-24) 0.0502 1.0515 0.1355 [0.806, 1.371] 0.7110
Age 45-54 (ref: 16-24) 0.4893 1.6311 0.1313 [1.261, 2.110] 0.0002 ***
Age 55-64 (ref: 16-24) 0.7224 2.0593 0.1286 [1.601, 2.649] 0.0000 ***
Age 65-74 (ref: 16-24) 0.9628 2.6190 0.1276 [2.040, 3.363] 0.0000 ***
Age 75+ (ref: 16-24) 1.0964 2.9933 0.1297 [2.322, 3.859] 0.0000 ***
Ethnicity: Black (ref: White) -0.2378 0.7883 0.1479 [0.590, 1.053] 0.1078
Ethnicity: Asian (ref: White) -0.3023 0.7391 0.1024 [0.605, 0.903] 0.0032 **
Ethnicity: Mixed (ref: White) -0.1601 0.8520 0.1641 [0.618, 1.175] 0.3293
Ethnicity: Other (ref: White) -0.4210 0.6564 0.2917 [0.371, 1.163] 0.1490
N=4,735 Log-lik=-5748.746 AIC=11525.492 BIC=11615.971 Alpha=0.3884
=== Full Coefficients: Model 3: + Full controls ===
Variable Coef IRR SE 95% CI p
----------------------------------------------------------------------------------------------------
GHQ-12 score 0.0921 1.0965 0.0055 [1.085, 1.108] 0.0000 ***
Female (ref: Male) 0.0154 1.0155 0.0377 [0.943, 1.093] 0.6824
Age 25-34 (ref: 16-24) -0.0368 0.9638 0.1387 [0.734, 1.265] 0.7906
Age 35-44 (ref: 16-24) 0.0507 1.0520 0.1351 [0.807, 1.371] 0.7077
Age 45-54 (ref: 16-24) 0.4547 1.5756 0.1309 [1.219, 2.037] 0.0005 ***
Age 55-64 (ref: 16-24) 0.6860 1.9857 0.1281 [1.545, 2.552] 0.0000 ***
Age 65-74 (ref: 16-24) 0.8958 2.4492 0.1284 [1.904, 3.150] 0.0000 ***
Age 75+ (ref: 16-24) 1.0066 2.7362 0.1326 [2.110, 3.548] 0.0000 ***
Ethnicity: Black (ref: White) -0.2758 0.7590 0.1479 [0.568, 1.014] 0.0622
Ethnicity: Asian (ref: White) -0.3069 0.7357 0.1049 [0.599, 0.904] 0.0034 **
Ethnicity: Mixed (ref: White) -0.1664 0.8467 0.1615 [0.617, 1.162] 0.3027
Ethnicity: Other (ref: White) -0.4417 0.6429 0.2893 [0.365, 1.134] 0.1268
Income Q2 (ref: Q1 lowest) 0.1112 1.1176 0.0557 [1.002, 1.246] 0.0458 *
Income Q3 (ref: Q1 lowest) -0.0119 0.9882 0.0608 [0.877, 1.113] 0.8455
Income Q4 (ref: Q1 lowest) -0.1217 0.8854 0.0647 [0.780, 1.005] 0.0598
Income Q5 (ref: Q1 lowest) -0.1595 0.8526 0.0713 [0.741, 0.980] 0.0252 *
NS-SEC: Intermediate (ref: Managerial) 0.0988 1.1039 0.0510 [0.999, 1.220] 0.0527
NS-SEC: Routine/manual (ref: Managerial) 0.1361 1.1458 0.0513 [1.036, 1.267] 0.0081 **
Education: HE below degree (ref: Degree) -0.0233 0.9770 0.0608 [0.867, 1.101] 0.7020
Education: A Level (ref: Degree) -0.1145 0.8918 0.0617 [0.790, 1.007] 0.0636
Education: O Level (ref: Degree) -0.0683 0.9339 0.0584 [0.833, 1.047] 0.2417
Education: NVQ1/CSE (ref: Degree) 0.0013 1.0013 0.1148 [0.800, 1.254] 0.9906
Education: Foreign/other (ref: Degree) -0.1818 0.8338 0.1880 [0.577, 1.205] 0.3334
Education: No qual (ref: Degree) -0.1256 0.8820 0.0668 [0.774, 1.005] 0.0601
Smoking: Ex-occasional (ref: Never) 0.0764 1.0794 0.0700 [0.941, 1.238] 0.2749
Smoking: Ex-regular (ref: Never) 0.1967 1.2173 0.0427 [1.120, 1.324] 0.0000 ***
Smoking: Current (ref: Never) 0.2175 1.2430 0.0600 [1.105, 1.398] 0.0003 ***
Alcohol: Non-drinker (ref: Lower risk) 0.1763 1.1928 0.0510 [1.079, 1.318] 0.0005 ***
Alcohol: Increased risk (ref: Lower risk) -0.1820 0.8336 0.0498 [0.756, 0.919] 0.0003 ***
Alcohol: Higher risk (ref: Lower risk) 0.0854 1.0892 0.0780 [0.935, 1.269] 0.2734
Region: North East (ref: London) 0.0168 1.0170 0.0873 [0.857, 1.207] 0.8473
Region: North West (ref: London) 0.1474 1.1588 0.0852 [0.980, 1.370] 0.0838
Region: Yorkshire (ref: London) -0.0342 0.9664 0.0879 [0.813, 1.148] 0.6973
Region: East Midlands (ref: London) 0.1169 1.1240 0.0908 [0.941, 1.343] 0.1982
Region: West Midlands (ref: London) -0.1211 0.8860 0.0939 [0.737, 1.065] 0.1973
Region: East of England (ref: London) 0.2382 1.2690 0.0842 [1.076, 1.497] 0.0046 **
Region: South East (ref: London) -0.0053 0.9948 0.0814 [0.848, 1.167] 0.9485
Region: South West (ref: London) 0.0141 1.0142 0.0855 [0.858, 1.199] 0.8690
IMD Q2 (ref: Q1 least deprived) 0.1384 1.1484 0.0550 [1.031, 1.279] 0.0119 *
IMD Q3 (ref: Q1 least deprived) 0.0650 1.0671 0.0587 [0.951, 1.197] 0.2684
IMD Q4 (ref: Q1 least deprived) 0.1342 1.1436 0.0603 [1.016, 1.287] 0.0261 *
IMD Q5 (ref: Q1 least deprived) 0.2253 1.2526 0.0638 [1.105, 1.419] 0.0004 ***
N=4,735 Log-lik=-5658.049 AIC=11404.098 BIC=11688.459 Alpha=0.3041
=== Full Coefficients: Model 4: + GHQ x IMD interaction ===
Variable Coef IRR SE 95% CI p
----------------------------------------------------------------------------------------------------
GHQ-12 score 0.0936 1.0981 0.0132 [1.070, 1.127] 0.0000 ***
Female (ref: Male) 0.0159 1.0160 0.0377 [0.944, 1.094] 0.6734
Age 25-34 (ref: 16-24) -0.0314 0.9691 0.1389 [0.738, 1.272] 0.8213
Age 35-44 (ref: 16-24) 0.0554 1.0570 0.1352 [0.811, 1.378] 0.6821
Age 45-54 (ref: 16-24) 0.4607 1.5851 0.1311 [1.226, 2.049] 0.0004 ***
Age 55-64 (ref: 16-24) 0.6927 1.9991 0.1284 [1.554, 2.571] 0.0000 ***
Age 65-74 (ref: 16-24) 0.9005 2.4609 0.1286 [1.913, 3.166] 0.0000 ***
Age 75+ (ref: 16-24) 1.0097 2.7448 0.1327 [2.116, 3.560] 0.0000 ***
Ethnicity: Black (ref: White) -0.2743 0.7601 0.1478 [0.569, 1.016] 0.0635
Ethnicity: Asian (ref: White) -0.3103 0.7333 0.1050 [0.597, 0.901] 0.0031 **
Ethnicity: Mixed (ref: White) -0.1694 0.8441 0.1615 [0.615, 1.158] 0.2941
Ethnicity: Other (ref: White) -0.4399 0.6441 0.2892 [0.365, 1.135] 0.1282
Income Q2 (ref: Q1 lowest) 0.1070 1.1129 0.0558 [0.998, 1.242] 0.0552
Income Q3 (ref: Q1 lowest) -0.0148 0.9853 0.0609 [0.874, 1.110] 0.8079
Income Q4 (ref: Q1 lowest) -0.1246 0.8828 0.0647 [0.778, 1.002] 0.0542
Income Q5 (ref: Q1 lowest) -0.1629 0.8497 0.0713 [0.739, 0.977] 0.0224 *
NS-SEC: Intermediate (ref: Managerial) 0.0986 1.1036 0.0511 [0.998, 1.220] 0.0536
NS-SEC: Routine/manual (ref: Managerial) 0.1351 1.1447 0.0514 [1.035, 1.266] 0.0086 **
Education: HE below degree (ref: Degree) -0.0239 0.9763 0.0609 [0.867, 1.100] 0.6941
Education: A Level (ref: Degree) -0.1150 0.8913 0.0618 [0.790, 1.006] 0.0628
Education: O Level (ref: Degree) -0.0686 0.9337 0.0584 [0.833, 1.047] 0.2403
Education: NVQ1/CSE (ref: Degree) 0.0024 1.0024 0.1147 [0.801, 1.255] 0.9834
Education: Foreign/other (ref: Degree) -0.1765 0.8382 0.1880 [0.580, 1.212] 0.3480
Education: No qual (ref: Degree) -0.1256 0.8820 0.0669 [0.774, 1.006] 0.0604
Smoking: Ex-occasional (ref: Never) 0.0783 1.0814 0.0701 [0.943, 1.241] 0.2641
Smoking: Ex-regular (ref: Never) 0.1982 1.2192 0.0427 [1.121, 1.326] 0.0000 ***
Smoking: Current (ref: Never) 0.2179 1.2434 0.0601 [1.105, 1.399] 0.0003 ***
Alcohol: Non-drinker (ref: Lower risk) 0.1781 1.1950 0.0511 [1.081, 1.321] 0.0005 ***
Alcohol: Increased risk (ref: Lower risk) -0.1817 0.8338 0.0498 [0.756, 0.919] 0.0003 ***
Alcohol: Higher risk (ref: Lower risk) 0.0875 1.0914 0.0780 [0.937, 1.272] 0.2620
Region: North East (ref: London) 0.0187 1.0189 0.0873 [0.859, 1.209] 0.8304
Region: North West (ref: London) 0.1485 1.1600 0.0853 [0.981, 1.371] 0.0817
Region: Yorkshire (ref: London) -0.0346 0.9660 0.0879 [0.813, 1.148] 0.6938
Region: East Midlands (ref: London) 0.1150 1.1219 0.0909 [0.939, 1.341] 0.2060
Region: West Midlands (ref: London) -0.1199 0.8870 0.0939 [0.738, 1.066] 0.2017
Region: East of England (ref: London) 0.2401 1.2714 0.0842 [1.078, 1.500] 0.0044 **
Region: South East (ref: London) -0.0049 0.9951 0.0814 [0.848, 1.167] 0.9518
Region: South West (ref: London) 0.0146 1.0147 0.0856 [0.858, 1.200] 0.8643
IMD Q2 (ref: Q1 least deprived) 0.1444 1.1553 0.0646 [1.018, 1.311] 0.0255 *
IMD Q3 (ref: Q1 least deprived) 0.0443 1.0453 0.0691 [0.913, 1.197] 0.5216
IMD Q4 (ref: Q1 least deprived) 0.1387 1.1488 0.0710 [1.000, 1.320] 0.0507
IMD Q5 (ref: Q1 least deprived) 0.2500 1.2840 0.0749 [1.109, 1.487] 0.0008 ***
GHQ x IMD Q2 (ref: Q1) -0.0032 0.9968 0.0177 [0.963, 1.032] 0.8552
GHQ x IMD Q3 (ref: Q1) 0.0086 1.0086 0.0178 [0.974, 1.044] 0.6301
GHQ x IMD Q4 (ref: Q1) -0.0024 0.9976 0.0179 [0.963, 1.033] 0.8941
GHQ x IMD Q5 (ref: Q1) -0.0090 0.9910 0.0172 [0.958, 1.025] 0.6009
N=4,735 Log-lik=-5657.445 AIC=11410.889 BIC=11721.101 Alpha=0.3035
Appendix D.2: Interaction Terms and Joint Wald Test¶
In [68]:
# ================================================================
# Appendix D.2: Interaction Terms and Tests of Effect Modification
# ================================================================
# --- Table D.2a: GHQ x IMD Interaction Coefficients ---
print("=== Table D.2a: GHQ x IMD Interaction Coefficients (Model 4) ===")
print("Note: ghq12scr = GHQ-12 slope in IMD Q1 (reference group).")
print(" Interaction terms = additional GHQ slope in each quintile vs Q1.")
print(f"\n{'Variable':<40} {'Coef':>8} {'IRR':>8} {'SE':>8} {'95% CI':>16} {'p':>8} {'':>5}")
print("-" * 90)
interact_vars = [
('GHQ-12 score (ref: IMD Q1)', 'ghq12scr'),
('GHQ x IMD Q2 (vs Q1)', 'ghq_x_imd2'),
('GHQ x IMD Q3 (vs Q1)', 'ghq_x_imd3'),
('GHQ x IMD Q4 (vs Q1)', 'ghq_x_imd4'),
('GHQ x IMD Q5 (vs Q1)', 'ghq_x_imd5'),
]
for label, var in interact_vars:
c = m4.params[var]
se = m4.bse[var]
p = m4.pvalues[var]
irr = np.exp(c)
ci_lo = np.exp(c - 1.96*se)
ci_hi = np.exp(c + 1.96*se)
stars = '***' if p<0.001 else '**' if p<0.01 else '*' if p<0.05 else ''
print(f"{label:<40} {c:>8.4f} {irr:>8.4f} {se:>8.4f} "
f"[{ci_lo:.3f}, {ci_hi:.3f}] {p:>8.4f} {stars:>5}")
# --- Table D.2b: Tests of Joint Significance ---
print("\n=== Table D.2b: Tests of Joint Significance of Interaction Terms ===")
print("H0: All GHQ x IMD interaction terms jointly equal zero")
print(" (i.e. the GHQ-12 effect does not vary across deprivation quintiles)")
# Likelihood Ratio Test
lr_stat = 2 * (m4.llf - m3.llf)
df_diff = len(inter)
lr_p = chi2_dist.sf(lr_stat, df_diff)
print(f"\nLikelihood Ratio Test (Model 4 vs Model 3):")
print(f" LR statistic: {lr_stat:.4f}")
print(f" Degrees of freedom:{df_diff}")
print(f" p-value: {lr_p:.4f}")
# Wald Test
inter_indices = [list(m4.params.index).index(v) for v in inter]
cov = np.array([[m4.cov_params()[i, j]
for j in inter_indices]
for i in inter_indices])
coefs = m4.params[inter].values
wald_stat = coefs @ np.linalg.inv(cov) @ coefs
wald_p = chi2_dist.sf(wald_stat, df_diff)
print(f"\nWald Test:")
print(f" Wald statistic: {wald_stat:.4f}")
print(f" Degrees of freedom:{df_diff}")
print(f" p-value: {wald_p:.4f}")
# Model comparison
print(f"\nModel Comparison (AIC):")
print(f" Model 3 (no interaction): AIC = {m3.aic:.3f}")
print(f" Model 4 (with interaction): AIC = {m4.aic:.3f}")
print(f" Delta AIC (M4 - M3): {m4.aic - m3.aic:.3f}")
print("\nNote: Interpret LR and Wald test results in conjunction.")
print(" See Appendix D.2 narrative for discussion.")
=== Table D.2a: GHQ x IMD Interaction Coefficients (Model 4) ===
Note: ghq12scr = GHQ-12 slope in IMD Q1 (reference group).
Interaction terms = additional GHQ slope in each quintile vs Q1.
Variable Coef IRR SE 95% CI p
------------------------------------------------------------------------------------------
GHQ-12 score (ref: IMD Q1) 0.0936 1.0981 0.0132 [1.070, 1.127] 0.0000 ***
GHQ x IMD Q2 (vs Q1) -0.0032 0.9968 0.0177 [0.963, 1.032] 0.8552
GHQ x IMD Q3 (vs Q1) 0.0086 1.0086 0.0178 [0.974, 1.044] 0.6301
GHQ x IMD Q4 (vs Q1) -0.0024 0.9976 0.0179 [0.963, 1.033] 0.8941
GHQ x IMD Q5 (vs Q1) -0.0090 0.9910 0.0172 [0.958, 1.025] 0.6009
=== Table D.2b: Tests of Joint Significance of Interaction Terms ===
H0: All GHQ x IMD interaction terms jointly equal zero
(i.e. the GHQ-12 effect does not vary across deprivation quintiles)
Likelihood Ratio Test (Model 4 vs Model 3):
LR statistic: 1.2090
Degrees of freedom:4
p-value: 0.8766
Wald Test:
Wald statistic: 1.2085
Degrees of freedom:4
p-value: 0.8767
Model Comparison (AIC):
Model 3 (no interaction): AIC = 11404.098
Model 4 (with interaction): AIC = 11410.889
Delta AIC (M4 - M3): 6.791
Note: Interpret LR and Wald test results in conjunction.
See Appendix D.2 narrative for discussion.
Appendix E: Marginal Effects¶
Appendix E.1: Predicted Condition Count Table¶
In [73]:
# ================================================================
# Appendix E.1: Predicted Condition Count Table
# ================================================================
coef_vars_m3 = [v for v in m3.params.index if v != 'alpha']
def avg_marginal_pred(ghq_val, imd_q, model):
r_temp = r.copy()
r_temp['ghq12scr'] = float(ghq_val)
for i in range(2, 6):
r_temp[f'imd_{i}'] = (1 if i == imd_q else 0)
X_temp = np.column_stack([
np.ones(len(r_temp)) if v == 'const'
else r_temp[v].values
for v in coef_vars_m3
])
coefs = model.params[coef_vars_m3].values
lp = X_temp @ coefs
return np.exp(lp).mean()
ghq_values = [0, 2, 4, 6, 8, 10, 12]
pred_matrix = {}
for ghq_val in ghq_values:
row = {}
for q in range(1, 6):
row[q] = avg_marginal_pred(ghq_val, q, m3)
pred_matrix[ghq_val] = row
# Table E.1a: Full — Q1, Q3, Q5 only
print("=== Table E.1a: Average Marginal Predicted Condition Count ===")
print("Model 3. All other covariates at individual observed values.")
print(f"\n{'GHQ Score':<12} {'IMD Q1':>12} {'IMD Q3':>12} {'IMD Q5':>12}")
print("-" * 50)
for ghq_val in ghq_values:
q1 = pred_matrix[ghq_val][1]
q3 = pred_matrix[ghq_val][3]
q5 = pred_matrix[ghq_val][5]
print(f"{ghq_val:<12} {q1:>12.3f} {q3:>12.3f} {q5:>12.3f}")
print(f"\n{'Absolute gap (GHQ=12 vs GHQ=0)':<35}", end='')
for q in [1, 3, 5]:
gap = pred_matrix[12][q] - pred_matrix[0][q]
print(f"{gap:>12.3f}", end='')
print()
# Table E.1b: Condensed — GHQ 0, 4, 8, 12 only
print("\n=== Table E.1b: Condensed Predicted Condition Count ===")
print("GHQ scores 0, 4, 8, 12 and IMD Q1, Q3, Q5 only.")
print(f"\n{'GHQ Score':<12} {'IMD Q1':>12} {'IMD Q3':>12} {'IMD Q5':>12}")
print("-" * 50)
for ghq_val in [0, 4, 8, 12]:
q1 = pred_matrix[ghq_val][1]
q3 = pred_matrix[ghq_val][3]
q5 = pred_matrix[ghq_val][5]
print(f"{ghq_val:<12} {q1:>12.3f} {q3:>12.3f} {q5:>12.3f}")
=== Table E.1a: Average Marginal Predicted Condition Count === Model 3. All other covariates at individual observed values. GHQ Score IMD Q1 IMD Q3 IMD Q5 -------------------------------------------------- 0 0.643 0.686 0.805 2 0.773 0.824 0.968 4 0.929 0.991 1.164 6 1.117 1.192 1.399 8 1.343 1.433 1.682 10 1.615 1.723 2.023 12 1.941 2.072 2.432 Absolute gap (GHQ=12 vs GHQ=0) 1.299 1.386 1.627 === Table E.1b: Condensed Predicted Condition Count === GHQ scores 0, 4, 8, 12 and IMD Q1, Q3, Q5 only. GHQ Score IMD Q1 IMD Q3 IMD Q5 -------------------------------------------------- 0 0.643 0.686 0.805 4 0.929 0.991 1.164 8 1.343 1.433 1.682 12 1.941 2.072 2.432
Appendix E.2: Marginal Effects Figure¶
In [74]:
# ================================================================
# Appendix E.2: Marginal Effects Figure
# ================================================================
ghq_fine = list(range(0, 13))
plot_imd = [1, 3, 5]
imd_colors = ['#1a9641', '#fdae61', '#d7191c']
imd_labels = ['Q1 (Least deprived)', 'Q3 (Middle)', 'Q5 (Most deprived)']
fig, ax = plt.subplots(figsize=(10, 6))
for q, col, lab in zip(plot_imd, imd_colors, imd_labels):
preds = [avg_marginal_pred(g, q, m3) for g in ghq_fine]
ax.plot(ghq_fine, preds,
color=col, linewidth=2.5,
marker='o', markersize=5,
label=lab)
ax.axvline(x=4, color='#999999', linestyle=':', linewidth=1.2, alpha=0.6)
ax.text(4.15, 0.75, 'Caseness\nthreshold', fontsize=8.5, color='#999999')
ax.set_xlabel('GHQ-12 Psychological Distress Score\n(0 = no distress, 12 = severe distress)',
fontsize=11, labelpad=8)
ax.set_ylabel('Average Predicted Physical Condition Count',
fontsize=11, labelpad=8)
ax.set_title('Figure E.2: Adjusted Physical Condition Count\n'
'by GHQ-12 Score and Area Deprivation, HSE 2022',
fontsize=12, fontweight='bold', pad=12)
ax.set_xticks(ghq_fine)
ax.legend(title='IMD Deprivation Quintile', fontsize=10,
title_fontsize=10, loc='upper left', framealpha=0.9)
ax.grid(alpha=0.25, linestyle='--')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
source_note_e2 = (
"Source: NHS England (2022) Health Survey for England, 2022 [data collection]. "
"UK Data Service, SN 9469. Average marginal predictions from Model 3 "
"(negative binomial regression, unweighted, N=4,735). "
"All other covariates retained at observed individual values."
)
wrapped_e2 = "\n".join(textwrap.wrap(source_note_e2, width=110))
fig.text(0.12, 0.01, wrapped_e2, fontsize=7.5, color="#555555")
plt.tight_layout(rect=[0, 0.10, 1, 1])
plt.savefig('figure_e2_marginal_effects.png', dpi=300, bbox_inches='tight')
plt.show()
print("Figure E.2 saved.")
Figure E.2 saved.
In [75]:
for bmi_v in ['bmival', 'bmivg3', 'bmisrg3']:
df[bmi_v] = clean(df[bmi_v])
valid = df[bmi_v].notna().sum()
pct = valid/len(df)*100
joint = df.dropna(subset=reg_vars+[bmi_v]).shape[0]
print(f'{bmi_v}: valid={valid:,} ({pct:.1f}%), joint sample={joint:,}')
bmival: valid=4,513 (58.4%), joint sample=3,293 bmivg3: valid=4,512 (58.4%), joint sample=3,293 bmisrg3: valid=2,469 (31.9%), joint sample=1,075
Appendix F: Robustness Checks¶
In [81]:
# ================================================================
# Appendix F: Robustness Checks
# ================================================================
from statsmodels.discrete.discrete_model import Poisson
from statsmodels.regression.linear_model import OLS
from statsmodels.discrete.discrete_model import Logit
from scipy.stats import chi2 as chi2_dist
# ================================================================
# F.1: Poisson vs Negative Binomial — LR Test
# ================================================================
print("=" * 70)
print("F.1: Poisson vs Negative Binomial — Likelihood Ratio Test")
print("=" * 70)
# Fit Poisson Model 3
def fit_poisson(xvars):
X = pd.DataFrame(
sm.add_constant(r[xvars].astype(float).values, has_constant='add'),
columns=['const'] + xvars
)
result = Poisson(r['condlcnt'].values, X.values).fit(
method='bfgs', maxiter=500, disp=False
)
param_names = ['const'] + xvars + []
result.params = pd.Series(result.params, index=['const'] + xvars)
result.bse = pd.Series(result.bse, index=['const'] + xvars)
result.pvalues = pd.Series(result.pvalues, index=['const'] + xvars)
return result
vars_m3 = ['ghq12scr'] + demo + ses + life + geo + imd
p3 = fit_poisson(vars_m3)
print(f"\nPoisson Model 3:")
print(f" Log-likelihood: {p3.llf:.3f}")
print(f" AIC: {p3.aic:.3f}")
print(f" GHQ-12 IRR: {np.exp(p3.params['ghq12scr']):.4f}")
print(f" GHQ-12 p: {p3.pvalues['ghq12scr']:.4f}")
print(f"\nNegative Binomial Model 3:")
print(f" Log-likelihood: {m3.llf:.3f}")
print(f" AIC: {m3.aic:.3f}")
print(f" GHQ-12 IRR: {np.exp(m3.params['ghq12scr']):.4f}")
print(f" GHQ-12 p: {m3.pvalues['ghq12scr']:.4f}")
# LR test: NB vs Poisson
# NB nests Poisson when alpha=0
# LR stat = 2*(llf_NB - llf_Poisson), df=1
lr_stat_f1 = 2 * (m3.llf - p3.llf)
lr_p_f1 = chi2_dist.sf(lr_stat_f1, 1) / 2 # one-sided
print(f"\nLikelihood Ratio Test (NB vs Poisson):")
print(f" LR statistic: {lr_stat_f1:.4f}")
print(f" p-value (one-sided): {lr_p_f1:.6f}")
print(f" Alpha (NB dispersion): {m3.params['alpha']:.4f}")
print(f"\nConclusion: {'Reject Poisson' if lr_p_f1 < 0.05 else 'Cannot reject Poisson'} in favour of NB")
# ================================================================
# F.2: OLS Robustness Check
# ================================================================
print("\n" + "=" * 70)
print("F.2: OLS Robustness Check")
print("=" * 70)
X_ols = sm.add_constant(
r[vars_m3].astype(float).values,
has_constant='add'
)
ols_result = OLS(r['condlcnt'].values.astype(float), X_ols).fit()
# Extract GHQ coefficient
ghq_idx = ['const'] + vars_m3
ghq_pos = ghq_idx.index('ghq12scr')
ols_coef = ols_result.params[ghq_pos]
ols_se = ols_result.bse[ghq_pos]
ols_p = ols_result.pvalues[ghq_pos]
ols_ci = ols_result.conf_int()[ghq_pos]
print(f"\nOLS Model (outcome: condlcnt, same controls as Model 3):")
print(f" GHQ-12 coefficient: {ols_coef:.4f}")
print(f" SE: {ols_se:.4f}")
print(f" 95% CI: [{ols_ci[0]:.4f}, {ols_ci[1]:.4f}]")
print(f" p-value: {ols_p:.4f}")
print(f" R-squared: {ols_result.rsquared:.4f}")
print(f" N: {int(ols_result.nobs):,}")
print(f"\nComparison with NB Model 3:")
print(f" NB IRR: {np.exp(m3.params['ghq12scr']):.4f}")
print(f" OLS coefficient: {ols_coef:.4f}")
print(f" Interpretation: A one-unit increase in GHQ-12 is associated with")
print(f" {ols_coef:.4f} additional conditions (OLS) vs {np.exp(m3.params['ghq12scr']):.4f}x IRR (NB)")
# ================================================================
# F.3: Logistic Regression — Binary Multimorbidity
# ================================================================
print("\n" + "=" * 70)
print("F.3: Logistic Regression — Binary Multimorbidity Outcome")
print("=" * 70)
def fit_logit(xvars):
X = pd.DataFrame(
sm.add_constant(r[xvars].astype(float).values, has_constant='add'),
columns=['const'] + xvars
)
result = Logit(r['multimorbid'].values.astype(float), X.values).fit(
method='bfgs', maxiter=500, disp=False
)
result.params = pd.Series(result.params, index=['const'] + xvars)
result.bse = pd.Series(result.bse, index=['const'] + xvars)
result.pvalues = pd.Series(result.pvalues, index=['const'] + xvars)
return result
logit_m3 = fit_logit(vars_m3)
c_l = logit_m3.params['ghq12scr']
se_l = logit_m3.bse['ghq12scr']
p_l = logit_m3.pvalues['ghq12scr']
or_l = np.exp(c_l)
ci_lo_l = np.exp(c_l - 1.96*se_l)
ci_hi_l = np.exp(c_l + 1.96*se_l)
print(f"\nLogistic Regression (outcome: multimorbid = 1 if condlcnt >= 2):")
print(f" GHQ-12 coefficient: {c_l:.4f}")
print(f" Odds Ratio: {or_l:.4f}")
print(f" 95% CI: [{ci_lo_l:.4f}, {ci_hi_l:.4f}]")
print(f" p-value: {p_l:.4f}")
print(f" N: {int(logit_m3.nobs):,}")
print(f"\nComparison with NB Model 3:")
print(f" NB IRR (condlcnt): {np.exp(m3.params['ghq12scr']):.4f}")
print(f" Logit OR (multimorbid): {or_l:.4f}")
print(f" Both statistically significant: {'Yes' if p_l < 0.001 else 'No'}")
# ================================================================
# F.4: BMI Sensitivity Analysis (nurse-measured bmival, N=3,293)
# ================================================================
print("\n" + "=" * 70)
print("F.4: BMI Sensitivity Analysis (nurse-measured bmival)")
print("=" * 70)
# Clean bmival
df_analytic['bmival_clean'] = df_analytic['bmival'].copy()
df_analytic.loc[df_analytic['bmival_clean'] < 0, 'bmival_clean'] = np.nan
# BMI subsample
bmi_vars = vars_m3 + ['bmival_clean']
r_bmi = df_analytic.dropna(subset=
['ghq12scr','condlcnt','qimd19','Sex','ag16g10',
'origin2','eqv5','nssec3','topqual3',
'cigst1_19','totalwug2_22','GOR1','bmival_clean']
).copy()
# Rebuild dummies for BMI subsample
for col in demo + ses + life + geo + imd:
if col not in r_bmi.columns:
r_bmi[col] = r[col]
print(f"\nBMI subsample N: {len(r_bmi):,}")
print(f"Main sample N: {len(r):,}")
print(f"Sample reduction: {len(r)-len(r_bmi):,} ({(len(r)-len(r_bmi))/len(r)*100:.1f}%)")
# Fit NB without BMI on subsample
def fit_nb_sub(xvars, data):
X = pd.DataFrame(
sm.add_constant(data[xvars].astype(float).values, has_constant='add'),
columns=['const'] + xvars
)
result = NegativeBinomial(data['condlcnt'].values, X.values).fit(
method='bfgs', maxiter=500, disp=False
)
result.params = pd.Series(result.params, index=['const'] + xvars + ['alpha'])
result.bse = pd.Series(result.bse, index=['const'] + xvars + ['alpha'])
result.pvalues = pd.Series(result.pvalues, index=['const'] + xvars + ['alpha'])
return result
m3_bmi_sub = fit_nb_sub(vars_m3, r_bmi)
m3_bmi_full = fit_nb_sub(vars_m3 + ['bmival_clean'], r_bmi)
print(f"\nNB without BMI (subsample N={len(r_bmi):,}):")
print(f" GHQ-12 IRR: {np.exp(m3_bmi_sub.params['ghq12scr']):.4f}")
print(f" p-value: {m3_bmi_sub.pvalues['ghq12scr']:.4f}")
print(f"\nNB with BMI (subsample N={len(r_bmi):,}):")
print(f" GHQ-12 IRR: {np.exp(m3_bmi_full.params['ghq12scr']):.4f}")
print(f" p-value: {m3_bmi_full.pvalues['ghq12scr']:.4f}")
print(f" BMI IRR: {np.exp(m3_bmi_full.params['bmival_clean']):.4f}")
print(f" BMI p: {m3_bmi_full.pvalues['bmival_clean']:.4f}")
print(f"\nComparison with main model (N={len(r):,}):")
print(f" Main model GHQ-12 IRR: {np.exp(m3.params['ghq12scr']):.4f}")
print(f" Subsample without BMI IRR: {np.exp(m3_bmi_sub.params['ghq12scr']):.4f}")
print(f" Subsample with BMI IRR: {np.exp(m3_bmi_full.params['ghq12scr']):.4f}")
====================================================================== F.1: Poisson vs Negative Binomial — Likelihood Ratio Test ====================================================================== Poisson Model 3: Log-likelihood: -5719.387 AIC: 11524.775 GHQ-12 IRR: 1.0920 GHQ-12 p: 0.0000 Negative Binomial Model 3: Log-likelihood: -5658.049 AIC: 11404.098 GHQ-12 IRR: 1.0965 GHQ-12 p: 0.0000 Likelihood Ratio Test (NB vs Poisson): LR statistic: 122.6766 p-value (one-sided): 0.000000 Alpha (NB dispersion): 0.3041 Conclusion: Reject Poisson in favour of NB ====================================================================== F.2: OLS Robustness Check ====================================================================== OLS Model (outcome: condlcnt, same controls as Model 3): GHQ-12 coefficient: 0.1024 SE: 0.0055 95% CI: [0.0917, 0.1132] p-value: 0.0000 R-squared: 0.1837 N: 4,735 Comparison with NB Model 3: NB IRR: 1.0965 OLS coefficient: 0.1024 Interpretation: A one-unit increase in GHQ-12 is associated with 0.1024 additional conditions (OLS) vs 1.0965x IRR (NB) ====================================================================== F.3: Logistic Regression — Binary Multimorbidity Outcome ====================================================================== Logistic Regression (outcome: multimorbid = 1 if condlcnt >= 2): GHQ-12 coefficient: 0.1627 Odds Ratio: 1.1767 95% CI: [1.1491, 1.2050] p-value: 0.0000 N: 4,735 Comparison with NB Model 3: NB IRR (condlcnt): 1.0965 Logit OR (multimorbid): 1.1767 Both statistically significant: Yes ====================================================================== F.4: BMI Sensitivity Analysis (nurse-measured bmival) ====================================================================== BMI subsample N: 3,293 Main sample N: 4,735 Sample reduction: 1,442 (30.5%) NB without BMI (subsample N=3,293): GHQ-12 IRR: 1.0929 p-value: 0.0000 NB with BMI (subsample N=3,293): GHQ-12 IRR: 1.0867 p-value: 0.0000 BMI IRR: 1.0344 BMI p: 0.0000 Comparison with main model (N=4,735): Main model GHQ-12 IRR: 1.0965 Subsample without BMI IRR: 1.0929 Subsample with BMI IRR: 1.0867