tooluniverse-epidemiological-analysis

流行病学数据分析

完整的观察性流行病学工作流——从研究问题到可直接发表的报告。每一步都编写并运行 Python 代码。不要只描述你“会怎么做”——直接动手完成。

步骤 1：提出研究问题（PECO 框架）

在接触数据之前，先定义P人群、E暴露、C对照、O结局。

• 人群：研究对象是谁？（例如，20-79 岁的成年人、III 期及以上癌症患者、ICU 收治患者）
• 暴露：研究什么因素？（例如，营养素摄入、药物治疗、基因突变、环境污染物）
• 对照：与什么进行比较？（例如，最低三分位组、未暴露组、野生型、安慰剂）
• 结局：研究什么健康事件？（例如，疾病发病率、生存时间、生物标志物水平、死亡率）

检查研究设计：研究问题是否要求考虑时间顺序？

• 横断面研究：某一时间点的患病率和关联
• 纵向研究/队列研究：发病率、因果推断、时间关系
• 病例对照研究：罕见结局、比值比（嵌套于队列中）
• 临床试验：随机对照下的干预效果

如果问题暗示因果关系（“X 是否会导致 Y？”），但只有横断面数据可用，请明确说明这一局限，并使用关联性表述继续分析。

步骤 2：查找并评估数据

使用 ToolUniverse 发现数据集，并查找既往研究使用过的数据：

# Search for relevant datasets — use find_tools to discover what's available
find_tools("dataset search")
find_tools("your domain keywords")  # e.g., "cancer genomics", "clinical trial", "survey health"

# Search literature for study precedents — papers cite their data sources
execute_tool("PubMed_search_articles", {"query": "[exposure] [outcome] [study design]", "max_results": 5})
execute_tool("EuropePMC_search_articles", {"query": "[exposure] [outcome] cohort", "limit": 5})

评估数据集是否适用：是否包含暴露变量、结局变量和关键混杂因素（年龄、性别以及领域特定因素）？样本量是否充足？

功效分析（在决定使用某数据集之前运行）：

from scipy.stats import norm
import numpy as np

def sample_size_logistic(p0, OR, alpha=0.05, power=0.80):
    """Minimum N for logistic regression detecting OR at given power."""
    p1 = (p0 * OR) / (1 - p0 + p0 * OR)
    z_a, z_b = norm.ppf(1 - alpha/2), norm.ppf(power)
    n = ((z_a + z_b)**2 * (1/(p0*(1-p0)) + 1/(p1*(1-p1)))) / (np.log(OR))**2
    return int(np.ceil(n))

print(f"Need N={sample_size_logistic(0.10, 1.5)} for OR=1.5 with 10% baseline prevalence")

步骤 3：下载并整理数据

以编程方式下载数据。根据数据源的格式调整加载代码。

import pandas as pd
import requests, io

# Generic download helper — adapt URL and format to your source
def download_and_parse(url, fmt="csv"):
    r = requests.get(url, timeout=120)
    content = io.BytesIO(r.content)
    if fmt == "xpt":
        return pd.read_sas(content, format="xport")
    elif fmt == "csv":
        return pd.read_csv(content)
    elif fmt == "tsv":
        return pd.read_csv(content, sep="\t")
    elif fmt == "stata":
        return pd.read_stata(content)
    elif fmt == "json":
        return pd.read_json(content)
    else:
        return pd.read_csv(content)  # default fallback

# Load and merge multiple files on shared ID column
df1 = download_and_parse(url1, fmt="xpt")
df2 = download_and_parse(url2, fmt="xpt")
df = df1.merge(df2, on="id_col", how="inner")

# Filter population (inclusion/exclusion criteria)
df = df[(df['age'] >= 20) & (df['age'] < 80)]

# Handle missing data
missing_pct = df.isnull().mean() * 100
print("Missing % per variable:\n", missing_pct[missing_pct > 0].sort_values(ascending=False))
# Decision: complete case if <5% missing; multiple imputation if 5-20%; drop variable if >20%

# Variable coding (adapt to your data)
df['age_group'] = pd.cut(df['age'], bins=[20,40,60,80], labels=['20-39','40-59','60-79'])
df['outcome_binary'] = (df['outcome_continuous'] >= threshold).astype(int)

调查权重：部分调查（NHANES、BRFSS、MEPS）需要使用抽样权重才能进行有效推断。请查阅调查文档。对于加权回归，请使用 statsmodels.stats.weightstats 或 linearmodels。

国际/中低收入国家（LMIC）人口健康：对于美国以外的 PECO 问题，DHSProgram_search_indicators/DHSProgram_get_data 可覆盖许多国家的同类国家级调查指标（生育率、孕产妇/儿童死亡率、营养、免疫接种、HIV），并可按国家和调查年份筛选——当研究人群不在美国时，这是类似 NHANES/BRFSS 的数据源。 欧盟健康数据集发现：对于聚焦欧洲的 PECO 问题，src/tooluniverse/data/euhealth_tools.json 提供约 20 个按主题划分的数据集发现工具（euhealthinfo_search_cancer、_search_deaths、_search_causes_of_death、_search_births、_search_infectious_diseases、_search_covid_19、_search_healthcare_expenditure、_search_diabetes_epidemiology_registry、_search_hospital_in_patient_data、_search_population_health_survey、_search_key_indicators_registries_surveys、_search_surveillance_mortality_rates、_search_surveillance、_search_primary_care_workforce、_search_obesity、_search_vaccination、_search_mental_health、_search_disability、_search_alcohol_tobacco_psychoactive_use、_search_cancer_registry），以及 euhealthinfo_deepdive。**这些工具用于查找相关欧洲数据集／登记库位于何处——返回的是数据集元数据（标题、落地页、关键词、许可、空间／语言覆盖范围），而不是原始统计数据本身。**它们都采用相同的参数格式：limit（默认值为 25）、country（全名或 ISO-3166 代码）、language（ISO 639-1 代码或全名）、term_override（用自定义文本替换工具内置的主题种子查询）、method（keyword／embedding／hybrid，默认值为 hybrid）、alpha（混合模式的权重比例）、top_k。响应格式（已实时验证）：{"results": [...]} 或一个裸数组，每条命中结果包含 uuid、title、landing_page、license、keywords、themes、language、spatial、snippet——将某条命中结果的 uuid（或直接使用搜索主题名称）传入 euhealthinfo_deepdive({"uuids": [...]}) 或 {"topic": "euhealthinfo_search_cancer", ...})，即可解析每个数据集分类后的外链（下载页面、资源门户）。

首次使用前必须完成设置：截至撰写本文时，这些工具会返回 {"warning": "EUHealth datastore not found locally (euhealth.db missing)...", "results": []}——已实时验证，并非猜测——除非先在本地同步底层数据集索引：

export HF_TOKEN=YOUR_HF_TOKEN
tu-datastore sync-hf download --repo "agenticx/tooluniverse-datastores" --collection euhealth --overwrite

在将空结果视为“没有匹配的数据集”之前，请检查响应中是否包含此警告——若 results 数组为空且出现数据存储警告，说明工具尚未完成设置，而不是没有找到任何内容。同步完成后，可将这组工具视为上文中 NHANES／DHS 发现模式在步骤 2 中的欧洲对应方案：先在此处找到合适的登记库／数据集，再进入步骤 3 实际下载并解析。

COVID-19 病例／死亡／疫苗接种时间序列：对于疫情监测类 PECO 问题，src/tooluniverse/data/diseasesh_tools.json 和 src/tooluniverse/data/disease_sh_ext_tools.json（Disease.sh API，数据来源为 JHU CSSE）可提供全球／国家／历史病例、死亡和康复数量，以及疫苗覆盖率——无需身份验证或设置（已实时验证，与上面的 EUHealth 不同）。DiseaseSh_get_global_stats（无参数）-> 全球总计；DiseaseSh_get_country_stats {"country": ...} -> 某个国家的当前快照；DiseaseSh_get_historical {"country": ..., "lastdays": ...} -> 每日累计时间线；DiseaseSH_get_vaccine_coverage {"country": ..., "lastdays": ...} -> 疫苗接种剂次时间线。请注意一个真实存在的近似重复项（已实时验证）：对于相同国家和日期范围，DiseaseSh_get_historical（来自 diseasesh_tools.json）和 DiseaseSH_get_covid_historical（来自 disease_sh_ext_tools.json）会从同一个底层 API 返回格式相同的时间线——任选其一即可；尽管它们位于不同的工具文件中，但并非两个独立的数据源。country: "all" 不是国家级工具的有效值（已实时验证：DiseaseSh_get_country_stats {"country": "all"} 会报错 "Country not found or doesn't have any cases"）——如需全球汇总，请使用 DiseaseSh_get_global_stats。 全球健康指标（WHO GHO）：对于跨国家比较，或需要使用 WHO 跟踪指标（预期寿命、非传染性疾病患病率、免疫接种覆盖率、死亡率、空气污染等——2000 多项指标）的 PECO 问题，WHOGHO_search_indicators/WHOGHO_get_indicator_data/WHOGHO_list_dimension_values 可直接查询 WHO 全球健康观察站 OData API，无需身份验证（已实时验证）。参数注意事项：这些工具使用 OData 筛选语法，而非普通关键词字符串——WHOGHO_search_indicators {"filter": "contains(IndicatorName,'life expectancy')", "top": 5}（不能用会报错的 {"query": ...}），以及 WHOGHO_get_indicator_data {"indicator_code": "WHOSIS_000001", "filter": "SpatialDim eq 'USA' and TimeDim eq 2019", "top": 5}，将范围限定为某个国家/年份（省略 filter 会返回涵盖多个国家/年份的完整数据序列，例如预期寿命有 12,936 行——务必筛选）。WHOGHO_list_dimension_values {"dimension_code": "COUNTRY"} 可查找有效的 SpatialDim 代码；传入 dimension_code: "Dimension" 可列出所有 195 个可筛选维度。已验证的实际数值：美国女性出生时预期寿命，2019 年 = 81.0 岁（WHOSIS_000001）。

实时环境暴露（空气质量）：对于包含空气污染暴露组的 PECO 问题，WAQI_get_air_quality {"city": "..."} 会返回指定城市/监测站当前由监测站测得的空气质量指数（已实时验证：柏林 -> AQI 87，data.iaqi 中包含主要污染物及各污染物的分指数），无需身份验证。OpenMeteo_get_air_quality {"latitude": ..., "longitude": ...} 则会返回任意经纬度位置逐小时的 PM2.5/PM10/CO/NO2/O3/SO2 模型浓度（不只是指数），还包括短期预报时段；当研究需要连续浓度变量而非分类 AQI，或位置附近没有 WAQI 监测站时，这是更好的数据源。将它们用作上述健康结局数据源的环境暴露对应数据（例如按城市-日期将 AQI/PM2.5 与结局数据集合并）——两者返回的历史时间序列都不会超过几天，因此若要获取回顾性队列的暴露史，请使用 OpenMeteo_get_historical_climate 或专门的空气质量存档。

美国县/社区层面的社会决定因素：health_disparities_get_svi_info/health_disparities_get_county_rankings_info（CDC/ATSDR 社会脆弱性指数和县级健康排名）——**已实时验证：这些是指针/元数据工具，而非数据检索工具。**尽管接受 state/county 参数，响应内容仍是指向来源网站的固定 access_url，外加对网站所提供内容的说明（下载 URL、数据年份、地理层级）——JSON 中不会返回各县的 SVI 分数或排名数值。使用这些工具查找数据集的存放位置及其涵盖内容，然后从返回的 URL 下载实际 CSV 进行分析；这与上述 euhealth 工具一样，属于“发现数据源，而非获取数据”的模式。

REST API 数据：对于 GDC（TCGA）、ClinicalTrials.gov 或 OpenTargets 等数据源，请通过分页遍历 API：

all_records = []
offset = 0
while True:
    resp = requests.get(f"{api_url}?offset={offset}&limit=500", timeout=30)
    batch = resp.json().get("data", [])
    if not batch:
        break
    all_records.extend(batch)
    offset += len(batch)
df = pd.DataFrame(all_records)

第 4 步：描述性统计（表 1）

# Table 1: mean +/- SD for continuous, N(%) for categorical, by exposure group
continuous_vars = ['age', 'bmi']  # adapt to your variables
for var in continuous_vars:
    print(df.groupby('exposure_group')[var].agg(['mean', 'std', 'count']))

categorical_vars = ['sex', 'race']  # adapt to your variables
for var in categorical_vars:
    print(pd.crosstab(df['exposure_group'], df[var], normalize='index') * 100)

检查分布：df[var].skew()、scipy.stats.shapiro()，以及用于识别离群值的直方图。

第 5 步：回归分析

逐步调整策略（为混杂因素提供证据）：

import statsmodels.formula.api as smf
import numpy as np

# Model 1: Unadjusted
m1 = smf.logit('outcome ~ exposure', data=df).fit(disp=0)

# Model 2: + demographics
m2 = smf.logit('outcome ~ exposure + age + sex + race', data=df).fit(disp=0)

# Model 3: + clinical factors
m3 = smf.logit('outcome ~ exposure + age + sex + race + bmi + smoking + alcohol', data=df).fit(disp=0)

# Report ORs with 95% CI
for name, model in [('Unadjusted', m1), ('Demographics', m2), ('Fully adjusted', m3)]:
    or_val = np.exp(model.params['exposure'])
    ci = np.exp(model.conf_int().loc['exposure'])
    print(f"{name}: OR={or_val:.2f} (95% CI: {ci[0]:.2f}-{ci[1]:.2f}), p={model.pvalues['exposure']:.4f}")

根据结局类型选择模型：

• 连续型结局：smf.ols()
• 二元结局：smf.logit()
• 有序分类：使用 statsmodels 中的 OrderedModel
• 生存时间结局：使用 lifelines 中的 CoxPHFitter
• 计数数据：smf.poisson() 或 smf.negativebinomial()

假设检验：

from statsmodels.stats.outliers_influence import variance_inflation_factor

# Multicollinearity (VIF > 5 is concerning, > 10 is severe)
X = df[['age', 'bmi', 'exposure']].dropna()
for i, col in enumerate(X.columns):
    print(f"VIF {col}: {variance_inflation_factor(X.values, i):.1f}")

第 6 步：敏感性分析

# Stratified analysis (effect modification)
for stratum_var in ['sex', 'age_group', 'race']:
    print(f"\n--- Stratified by {stratum_var} ---")
    for level, sub in df.groupby(stratum_var):
        if len(sub) < 50: continue
        try:
            m = smf.logit('outcome ~ exposure + age + bmi', data=sub).fit(disp=0)
            or_val = np.exp(m.params['exposure'])
            ci = np.exp(m.conf_int().loc['exposure'])
            print(f"  {level}: OR={or_val:.2f} ({ci[0]:.2f}-{ci[1]:.2f}), p={m.pvalues['exposure']:.3f}, N={len(sub)}")
        except: print(f"  {level}: model failed (N={len(sub)})")

# Exclude outliers (+/- 3 SD) and re-run
df_no_outliers = df[np.abs(stats.zscore(df['exposure'].dropna())) < 3]
m_robust = smf.logit('outcome ~ exposure + age + sex + bmi', data=df_no_outliers).fit(disp=0)

# Confounder-adjusted exposure (residual method, e.g., energy-adjusted nutrient intake)
# Use when exposure correlates strongly with a confounder (total calories, body size, etc.)
adj_model = smf.ols('exposure ~ confounder', data=df).fit()
df['exposure_adj'] = adj_model.resid

# Multiple comparisons note
n_tests = 5  # number of exposure-outcome pairs tested
bonferroni_threshold = 0.05 / n_tests

步骤 7：生物学解读（ToolUniverse 的优势）

ToolUniverse 的价值不仅限于统计软件包。在发现统计关联后，还要探究其生物学合理性。使用 find_tools 来查找适合你所在领域的工具。

# 1. Literature evidence — search for mechanism connecting exposure to outcome
execute_tool("PubMed_search_articles", {"query": "[exposure] [outcome] mechanism", "max_results": 5})
execute_tool("EuropePMC_search_articles", {"query": "[exposure] [outcome] mouse model in vivo", "limit": 5})

# 2. Pathway/molecular context — discover tools for your exposure type
find_tools("[exposure type] pathway")       # e.g., "nutrient pathway", "drug target", "chemical toxicology"
find_tools("[outcome type] gene disease")   # e.g., "diabetes gene", "cancer survival", "cardiac risk"

# 3. Gene-disease evidence (if a gene/variant is involved)
find_tools("gene disease association")
find_tools("variant functional annotation")

# 4. Drug/chemical mechanisms (if a drug or chemical is the exposure)
find_tools("drug mechanism target")
find_tools("chemical gene interaction")

模式：暴露 X →（涉及什么分子通路？）→（涉及什么生物过程？）→ 结局 Y。使用 ToolUniverse 工具补全中间步骤。这样可以将统计关联转化为生物学上合理的假设。

步骤 8：可视化

建议绘制的关键图表（使用 matplotlib）：

• 森林图：分层 OR，包含 95% CI；在 OR=1 处添加垂直参考线，x 轴使用对数刻度
• 剂量-反应曲线：x 轴为暴露四分位数的中位数，y 轴为结局患病率/均值
• DAG：显示假定因果结构（暴露、混杂因素、结局）的有向无环图
• 散点图 + 回归线：针对连续结局，使用 sns.regplot()

步骤 9：报告结构

按以下顺序撰写最终报告：

1. 背景：已知信息，以及本分析要解决的研究空白（引用步骤 7 中的 PubMed 检索）
2. 方法：PECO、数据来源、纳入/排除标准、变量定义、统计方法、敏感性分析
3. 结果：表 1、未调整和调整后的 OR/HR、分层结果、剂量-反应、敏感性检查
4. 讨论：与既往文献比较（PubMed）、生物学合理性（ToolUniverse 的通路/机制发现）、效应量的临床意义
5. 局限性：研究设计的限制、未测量的混杂因素、选择偏倚、测量误差、普遍适用性

始终应说明的关键局限性：

• 横断面设计无法确定时间先后关系（如适用）
• 未测量变量导致的残余混杂
• 自我报告的暴露数据可能存在回忆偏倚
• 调查权重可能无法充分校正无应答偏倚

完整性检查清单

完成任何流行病学分析前，请确认：

• 已定义并记录 PECO
• 研究设计与研究问题相匹配
• 样本量充足（已进行效能分析）
• 已报告并处理缺失数据
• 已生成表 1
• 已报告未调整和调整后的模型
• 已论证混杂因素的选择依据（而非仅依据统计结果）
• 已检查假设（VIF、线性、模型拟合）
• 至少进行了 1 项敏感性分析
• 已通过 ToolUniverse 调查生物学合理性
• 已结合临床背景解读效应量（而不只是 p 值）
• 已明确说明局限性