workflows/clinical-trial-pipeline/SKILL.md
End-to-end clinical trial analysis workflow from CDISC SDTM/ADaM loading through ICH E9(R1) estimand-driven primary analysis to CONSORT 2025 regulatory-compliant reporting. Covers data preparation, FDA 2023 marginal vs conditional logistic regression, categorical tests with Boschloo, modern HTE/subgroup methods, missing-data sensitivity (MMRM, reference-based MI, Permutt tipping point), graphical multiplicity (Bretz-Maurer), survival analysis (Cox/RMST/competing risks) when applicable, and Table 1. Use when performing a complete analysis of clinical trial data.
npx skillsauth add GPTomics/bioSkills bio-workflows-clinical-trial-pipelineInstall this skill globally with one command. Works with Claude Code, Cursor, and Windsurf.
3 of 9 scanners reported clean
Some scanners were skipped, did not run, or reported a non-clean status. Review each row below.
Reference examples tested with: statsmodels 0.14+, scipy 1.12+, tableone 0.9+, pyreadstat 1.2+, pandas 2.1+, numpy 1.26+, matplotlib 3.8+
Before using code patterns, verify installed versions match. If versions differ:
pip show <package> then help(module.function) to check signaturesIf code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
"Analyze my clinical trial data end to end" -> Load CDISC domain tables, prepare a subject-level analysis dataset, run primary statistical models, perform subgroup analyses, and generate regulatory-compliant tables and figures.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step.
A clinical-trial analysis is a commit-then-execute pipeline: its trustworthiness is decided by whether the analysis was locked BEFORE the data was seen, and every seam failure is a wrong-but-silent handoff that answers a different question with no error thrown.
| Commitment | Consequence inherited downstream | |------------|----------------------------------| | The estimand (5 ICH E9(R1) attributes, SAP-locked before unblinding) | Which question the trial answers; the estimator must match its ICE strategy | | The estimator matched to the ICE strategy | Whether the analysis answers the locked question; a mismatch is silent | | Analysis population set (ITT/FAS primary, PP sensitivity, Safety as-treated) a priori | Causal validity; ITT is randomization-protected, PP and subgroups are not | | SDTM -> ADaM derivation (pre-specified, one-directional; ADTTE CNSR convention) | Every downstream model; CDISC CNSR=0 means EVENT, opposite of R/Python survival packages |
Before executing any analysis step, establish the causal framework. For an RCT, randomization justifies causal interpretation of the primary analysis, but subgroup analyses and observational comparisons within the trial (e.g., adherence effects) do not inherit this protection. Key decisions requiring scientific judgment at each step: (1) data preparation -- which aggregation strategy matches the estimand, (2) covariate selection -- include confounders and prognostic factors from the SAP, exclude mediators and colliders, (3) subgroup analysis -- test only biologically motivated interactions, (4) missing data -- link DS domain reasons to the assumed mechanism before choosing a method. The workflow below provides the technical steps; the scientific reasoning at each decision point determines whether the results are valid.
CDISC Domain Files (DM, AE, EX, LB)
|
v
[1. Data Preparation] ----> Subject-level dataset with outcomes and covariates
|
v
[2. Table 1] ------------> Baseline characteristics by treatment arm
|
v
[3. Primary Analysis] ---> Marginal RD via g-computation (conditional OR supportive)
|
v
[4. Categorical Tests] --> Chi-square / Fisher's exact for key associations
|
v
[5. Subgroup Analysis] --> Interaction terms, stratified ORs, forest plot
|
v
[6. Missing Data] -------> Multiple imputation sensitivity analysis
|
v
Results tables and figures
Goal: Create a single subject-level analysis dataset from CDISC domain tables.
Approach: Load domain files, aggregate event-level data to one row per subject, merge on USUBJID, and code the outcome variable.
import pandas as pd
import pyreadstat
dm, _ = pyreadstat.read_xport('dm.xpt')
ae, _ = pyreadstat.read_xport('ae.xpt')
# Aggregate: did each subject have the target adverse event?
target_ae = ae[ae['AEDECOD'] == 'COVID-19'].copy()
severity_map = {'MILD': 1, 'MODERATE': 2, 'SEVERE': 3, 'LIFE THREATENING': 4, 'FATAL': 5}
target_ae['AESEV_NUM'] = target_ae['AESEV'].map(severity_map)
had_event = target_ae.groupby('USUBJID')['AESEV_NUM'].max().reset_index()
had_event.columns = ['USUBJID', 'EVENT_SEVERITY']
analysis = dm[['USUBJID', 'ARM', 'ARMCD', 'AGE', 'SEX']].merge(had_event, on='USUBJID', how='left')
analysis['HAD_EVENT'] = analysis['EVENT_SEVERITY'].notna().astype(int)
analysis['TREATMENT'] = (analysis['ARMCD'] != 'PLACEBO').astype(int)
QC Checkpoint: Verify one row per USUBJID, no unexpected duplicates, treatment arms are present and reasonably balanced.
assert analysis['USUBJID'].is_unique, 'Duplicate subjects detected'
print(analysis['ARM'].value_counts())
Goal: Summarize demographics and baseline variables by treatment arm.
Approach: Use TableOne to generate a baseline table by arm with standardized mean differences and explicit missingness. Omit the baseline p-value column: in a randomized trial any imbalance is by definition due to chance, so a baseline p-value tests a null already known to be true (Senn 1994; CONSORT 2010/2025). Report SMD for balance instead. Table-construction and export mechanics (gtsummary/tableone, Word export, gene-symbol-safe supplements) live in reporting/publication-tables.
from tableone import TableOne
columns = ['AGE', 'SEX', 'RACE']
categorical = ['SEX', 'RACE']
table1 = TableOne(analysis, columns=columns, categorical=categorical,
groupby='ARM', pval=False, smd=True, missing=True)
print(table1.tabulate(tablefmt='github'))
Interpret SMD > 0.1 as meaningful imbalance. The response to a worrying imbalance on a prognostic covariate is to adjust for it (a pre-specified ANCOVA/model covariate), not to test it.
Goal: Estimate the treatment effect on the binary outcome as an adjusted odds ratio.
Approach: Fit a logistic regression with explicit reference category and clinically relevant covariates, then exponentiate coefficients to obtain ORs.
import statsmodels.formula.api as smf
import numpy as np
model = smf.logit(
'HAD_EVENT ~ C(ARM, Treatment(reference="Placebo")) + AGE + C(SEX)',
data=analysis
).fit()
or_table = pd.DataFrame({
'OR': np.exp(model.params),
'Lower_CI': np.exp(model.conf_int()[0]),
'Upper_CI': np.exp(model.conf_int()[1]),
'p_value': model.pvalues
})
print(or_table)
print(f'McFadden pseudo-R2: {model.prsquared:.4f}')
The logistic fit above yields the CONDITIONAL (adjusted) OR. When the estimand's summary measure is a marginal risk difference (the FDA 2023 primary for a binary endpoint), do NOT report this conditional OR as the primary effect. Compute the MARGINAL risk difference by g-computation: fit the covariate-adjusted model, predict each subject's outcome probability under both arms, average within arm, and difference; bootstrap the whole fit-and-predict for the CI. examples/clinical_trial_pipeline.py implements this; see clinical-biostatistics/logistic-regression for the estimator's assumptions and variance options. Report the conditional OR as supportive. The two differ by non-collapsibility and answer different questions.
QC Checkpoint: Verify model converged (no warnings), check for separation (coefficients > 10 or SE > 100), report pseudo-R-squared (McFadden > 0.2 is excellent; do not compare across pseudo-R2 types). Confirm the reported PRIMARY effect matches the estimand's summary measure (marginal RD via g-computation for a marginal estimand), not whichever the model emits by default.
Goal: Test the crude association between treatment and outcome using contingency tables.
Approach: Build a 2x2 table, check expected cell counts, and choose chi-square or Fisher's exact accordingly.
from scipy.stats import chi2_contingency, fisher_exact
table = pd.crosstab(analysis['ARM'], analysis['HAD_EVENT'])
chi2, p, dof, expected = chi2_contingency(table, correction=False)
if (expected < 5).any():
_, p = fisher_exact(table.values)
print(f'Fisher exact p = {p:.4f}')
else:
print(f'Chi-square p = {p:.4f} (chi2 = {chi2:.2f}, dof = {dof})')
Goal: Test whether the treatment effect varies across pre-specified subgroups.
Approach: Fit a model with an interaction term and test it with a single interaction LR test. Subgroup-specific ORs are DESCRIPTIVE ONLY -- do not correct their individual p-values, do not interpret them; the alpha allocated to the subgroup family is handled by the pre-specified gMCP graph, not by a post-hoc correction. Visualize with a forest plot.
import matplotlib.pyplot as plt
# The HTE test is the treatment-by-subgroup INTERACTION, not per-subgroup significance.
# Fit the main-effects (restricted) and interaction (full) models, then LR-test the interaction.
main_model = smf.logit(
'HAD_EVENT ~ C(ARM, Treatment(reference="Placebo")) + C(SUBGROUP)',
data=analysis
).fit(disp=0)
interaction_model = smf.logit(
'HAD_EVENT ~ C(ARM, Treatment(reference="Placebo")) * C(SUBGROUP)',
data=analysis
).fit(disp=0)
# compare_lr_test exists only on linear-model results, NOT LogitResults -- compute the LR test by hand.
from scipy.stats import chi2
lr_stat = 2 * (interaction_model.llf - main_model.llf)
df_diff = int(interaction_model.df_model - main_model.df_model)
interaction_pval = chi2.sf(lr_stat, df_diff)
print(f'Treatment-by-subgroup interaction: LR chi2={lr_stat:.2f}, df={df_diff}, p={interaction_pval:.3f}')
# Subgroup-specific ORs are DESCRIPTIVE ONLY (to draw the forest plot). Do NOT interpret their
# individual p-values as evidence of a subgroup effect -- that is the invalid per-subgroup
# significance pattern; the interaction p-value above is the only valid HTE test.
labels, ors, lowers, uppers = [], [], [], []
for group in analysis['SUBGROUP'].unique():
sub = analysis[analysis['SUBGROUP'] == group]
sub_model = smf.logit(
'HAD_EVENT ~ C(ARM, Treatment(reference="Placebo"))',
data=sub
).fit(disp=0)
or_val = np.exp(sub_model.params.iloc[1])
ci = np.exp(sub_model.conf_int().iloc[1])
labels.append(group)
ors.append(or_val)
lowers.append(ci[0])
uppers.append(ci[1])
# Forest plot
fig, ax = plt.subplots(figsize=(8, 5))
y_pos = range(len(labels))
ax.errorbar(ors, y_pos,
xerr=[np.array(ors) - np.array(lowers), np.array(uppers) - np.array(ors)],
fmt='D', color='black', capsize=3, markersize=5)
ax.axvline(x=1.0, color='gray', linestyle='--', linewidth=0.8)
ax.set_yticks(y_pos)
ax.set_yticklabels(labels)
ax.set_xlabel('Odds Ratio (95% CI)')
ax.set_xscale('log')
plt.tight_layout()
plt.savefig('forest_plot.png', dpi=150)
QC Checkpoint: Interaction p-value reported from a single LR test. Alpha for the subgroup family allocated in the pre-specified gMCP graph (not a post-hoc p-value correction on the descriptive subgroup ORs). Forest plot shows overall estimate for context.
Goal: Assess robustness of the primary result under the pre-specified ICE strategy with both MAR primary and MNAR sensitivity analyses.
Approach: First examine DS (Disposition) domain for differential dropout patterns; if dropout differs by arm, MAR is suspect and reference-based MI is required as primary. Otherwise, fit MMRM under MAR with Rubin's-rules pooling for continuous endpoints, or g-computation with bootstrap for binary. Always run Permutt 2016 tipping-point sensitivity in residual SD units.
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
n_imputations = 20 # practical starting count; von Hippel 2020 shows required m scales with the fraction of missing information (two-stage rule), so raise it when FMI is high
# Impute AGE jointly WITH its predictors (SEX, TREATMENT, HAD_EVENT). A single-column imputer has
# no predictors, so sample_posterior draws are identical -> between-imputation variance = 0 and the
# MI collapses to complete-case. Including correlated columns makes the posterior draws actually vary.
impute_cols = ['AGE', 'TREATMENT', 'HAD_EVENT'] # numeric only; SEX ('M'/'F') would break IterativeImputer and is restored below
mi_data = analysis.dropna(subset=['HAD_EVENT', 'TREATMENT']).copy()
results = []
for i in range(n_imputations):
imputer = IterativeImputer(max_iter=10, random_state=i, sample_posterior=True)
imputed_cov = pd.DataFrame(imputer.fit_transform(mi_data[impute_cols]),
columns=impute_cols, index=mi_data.index)
# Only AGE had missings; restore the observed discrete columns so they stay integer-valued.
imputed_cov['HAD_EVENT'] = mi_data['HAD_EVENT'].values
imputed_cov['TREATMENT'] = mi_data['TREATMENT'].values
imputed_cov['SEX'] = mi_data['SEX'].values
# Mirror the primary ADJUSTED model's RHS (AGE + SEX) so the pooled OR is comparable to it.
model_imp = smf.logit('HAD_EVENT ~ TREATMENT + AGE + C(SEX)', data=imputed_cov).fit(disp=0)
results.append({'coef': model_imp.params['TREATMENT'], 'se': model_imp.bse['TREATMENT']})
pooled_coef = np.mean([r['coef'] for r in results])
within_var = np.mean([r['se']**2 for r in results])
between_var = np.var([r['coef'] for r in results], ddof=1)
total_var = within_var + (1 + 1/n_imputations) * between_var
pooled_or = np.exp(pooled_coef)
pooled_ci = (np.exp(pooled_coef - 1.96 * np.sqrt(total_var)),
np.exp(pooled_coef + 1.96 * np.sqrt(total_var)))
print(f'Pooled OR: {pooled_or:.3f} ({pooled_ci[0]:.3f}-{pooled_ci[1]:.3f})')
QC Checkpoint: Compare pooled OR and CI with the complete-case primary analysis. Large discrepancies suggest missing data may not be MCAR. Document the comparison.
| Symptom | Cause | Fix | |---------|-------|-----| | The analysis answers a different question | Estimand-estimator mismatch (e.g. a hypothetical estimand analyzed treatment-policy) | Derive the estimator FROM the ICE strategy; MMRM/MI/g-computation are consequences of attribute (iv), not free choices | | Spurious subgroup effect | Subgroups/sensitivity chosen after seeing the data | Pre-specify all subgroups + the multiplicity graph in the SAP; post-hoc is hypothesis-generating only | | Effect loses causal validity | PP swapped in as primary because it "looks cleaner" | ITT/FAS primary (randomization-protected); PP is sensitivity; define both a priori | | Survival results inverted | ADTTE CNSR sign flip (CDISC CNSR=0 means EVENT) | Convert the censoring indicator before passing to lifelines/survival | | False "imbalance" conclusions | Baseline characteristics tested with p-values | Report SMD (>0.1 notable); adjust prognostic imbalance via pre-specified ANCOVA | | Primary effect is the wrong parameter | Conditional OR reported as the marginal effect (non-collapsibility) | Marginal risk difference via g-computation is primary (FDA 2023); conditional OR supportive | | Over-conservative inference (CIs too wide, power lost) | Stratification/randomization factors omitted from the analysis model | Include the stratification factors; omitting them biases SEs upward and Type-I below nominal (Kahan-Morris 2012) |
This pipeline covers the typical binary-endpoint RCT workflow. For specific designs, add the corresponding specialized skill:
tools
End-to-end CLIP-seq pipeline from FASTQ to ENCODE-compliant binding sites, single-nucleotide crosslink maps, annotation, motifs, and (optionally) differential binding. Use when running the full Yeo lab eCLIP / iCLIP / iCLIP2 / iCLIP3 / irCLIP / PAR-CLIP analysis with SMInput control, protocol-specific UMI extraction, ENCODE STAR parameters, CLIPper or Skipper peak calling with stringent log2 FC and -log10 p thresholds, IDR rescue and self-consistency QC, and downstream motif registration with mCross or PEKA.
development
Detect, date, and contextualize whole-genome duplication (WGD / paleopolyploidy) events using wgd v2 (Chen et al 2024), KsRates (Sensalari 2022 substitution-rate-corrected Ks dating), DupGen_finder (Qiao 2019), MAPS (Li 2018 phylogenomic), POInT (Conant 2008 ordered-block), SLEDGe (2024 ML-based), Whale.jl (Bayesian DL+WGD), and synteny-anchored paranome construction. Use when identifying ancient polyploidy from Ks distributions and synteny block analysis, positioning WGD events relative to speciation, distinguishing tandem from segmental from WGD duplications, dating the 2R/3R vertebrate / fish / salmonid WGDs, building paranome and Ks-age mixture models, applying KsRates substitution-rate correction across lineages, or testing alternative biased-fractionation / dosage-balance models post-WGD.
tools
Build whole-genome alignments using Progressive Cactus (Armstrong 2020 reference-free clade-level WGA), Minigraph-Cactus (Hickey 2024 pangenome-aware), LASTZ chain/net (UCSC pipeline), MUMmer4 (Marçais 2018 pairwise), minimap2 -x asm5/10/20 (Li 2018 fast pairwise), AnchorWave (Song 2022 WGD-aware), and Mauve / progressiveMauve (bacterial). Operates the HAL toolkit (Hickey 2013) for downstream extraction including halSynteny, halLiftover, halBranchMutations, and hal2maf. Use when constructing multi-species alignments for comparative-annotation projection (TOGA), synteny detection, conservation analyses (phyloP / PhastCons), or pangenome graph construction; selecting between reference-free (Cactus) and reference-anchored (LASTZ chains/nets) approaches; tuning sensitivity for closely vs distantly related genomes; or producing HAL files for genome-wide downstream tools.
development
Detect syntenic blocks and structural rearrangements between genomes using MCScanX (Wang 2012), JCVI/MCScan (Tang 2008 Python), GENESPACE (Lovell 2022) for orthology-anchored riparian visualization, SyRI for structural variation, AnchorWave for sequence-level synteny, i-ADHoRe 3.0 for highly diverged species, SynNet for synteny networks, and ntSynt for multi-genome macrosynteny. Use when identifying collinear gene blocks across species, distinguishing macrosynteny from microsynteny, detecting inversions/translocations/duplications, anchoring orthology in WGD lineages, producing publication riparian plots, computing synteny block age via Ks (cross-references whole-genome-duplication), or running synteny-aware ortholog inference in polyploids.